Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 19 for palm/trunk


Ignore:
Timestamp:
Feb 23, 2007 4:53:48 AM (18 years ago)
Author:
raasch
Message:

preliminary version of modified boundary conditions at top

Location:
palm/trunk/SOURCE
Files:
23 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/advec_s_pw.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt
    77!
    88! Former revisions:
    … …  
    5858       DO  i = nxl, nxr
    5959          DO  j = nys, nyn
    60              DO  k = nzb_s_inner(j,i)+1, nzt-1
     60             DO  k = nzb_s_inner(j,i)+1, nzt
    6161                tend(k,j,i) = tend(k,j,i)                                      &
    6262              -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
    … …  
    9393
    9494
    95        DO  k = nzb_s_inner(j,i)+1, nzt-1
     95       DO  k = nzb_s_inner(j,i)+1, nzt
    9696          tend(k,j,i) = tend(k,j,i)                                            &
    9797              -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &

  • palm/trunk/SOURCE/boundary_conds.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
     7! values are now calculated by the prognostic equation,
     8! Dirichlet and zero gradient condition for pt established at top boundary
    79!
    810! Former revisions:
    … …  
    100102!
    101103!--    Temperature at top boundary
    102        IF ( ibc_pt_t == 1 )  THEN
    103           pt(nzt,:,:)   = pt(nzt-1,:,:) + bc_pt_t_val * dzu(nzt)
     104       IF ( ibc_pt_t == 0 )  THEN
     105          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     106             pt(nzt+1,:,:) = pt_m(nzt+1,:,:)
     107          ELSE
     108             pt(nzt+1,:,:) = pt_p(nzt+1,:,:)  ! pt_m not used for Runge-Kutta
     109          ENDIF
     110       ELSEIF ( ibc_pt_t == 1 )  THEN
     111          pt(nzt+1,:,:) = pt(nzt,:,:)
     112       ELSEIF ( ibc_pt_t == 2 )  THEN
    104113          pt(nzt+1,:,:) = pt(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
    105114       ENDIF
    … …  
    114123             ENDDO
    115124          ENDDO
    116           e(nzt,:,:)   = e(nzt-1,:,:)
    117125          e(nzt+1,:,:) = e(nzt,:,:)
    118126       ENDIF
    … …  
    147155!
    148156!--       Top boundary
    149           q(nzt,:,:)   = q(nzt-1,:,:) + bc_q_t_val * dzu(nzt)
    150157          q(nzt+1,:,:) = q(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
    151158       ENDIF

  • palm/trunk/SOURCE/calc_liquid_water_content.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Old comment removed
    77!
    88! Former revisions:
    … …  
    7979    ENDDO
    8080   
    81 !
    82 !-- Setting boundary values of ql
    83 !    CALL exchange_horiz( ql, 0, 0 )   
    84 
    8581 END SUBROUTINE calc_liquid_water_content

  • palm/trunk/SOURCE/calc_precipitation.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt
    77!
    88! Former revisions:
    … …  
    5454       DO  i = nxl, nxr
    5555          DO  j = nys, nyn
    56              DO  k = nzb_2d(j,i)+1, nzt-1   
     56             DO  k = nzb_2d(j,i)+1, nzt
    5757
    5858                IF ( ql(k,j,i) > ql_crit )  THEN
    … …  
    8787
    8888
    89        DO  k = nzb_2d(j,i)+1, nzt-1   
     89       DO  k = nzb_2d(j,i)+1, nzt
    9090
    9191          IF ( ql(k,j,i) > ql_crit )  THEN

  • palm/trunk/SOURCE/check_parameters.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Temperature and humidity gradients at top are now calculated for nzt+1,
     7! top_heatflux and respective boundary condition bc_pt_t is checked
    78!
    89! Former revisions:
    … …  
    483484!--    Store temperature gradient at the top boundary for possile Neumann
    484485!--    boundary condition
    485        bc_pt_t_val = ( pt_init(nzt) - pt_init(nzt-1) ) / dzu(nzt)
     486       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
    486487
    487488!
    … …  
    539540!--       Store humidity gradient at the top boundary for possile Neumann
    540541!--       boundary condition
    541           bc_q_t_val = ( q_init(nzt) - q_init(nzt-1) ) / dzu(nzt)
     542          bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
    542543
    543544       ENDIF
    … …  
    795796    ELSEIF ( bc_pt_t == 'neumann' )  THEN
    796797       ibc_pt_t = 1
     798    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
     799       ibc_pt_t = 2
    797800    ELSE
    798801       IF ( myid == 0 )  THEN
    … …  
    804807
    805808    IF ( surface_heatflux == 9999999.9 )  constant_heatflux = .FALSE.
     809    IF ( top_heatflux     == 9999999.9 )  THEN
     810       constant_top_heatflux = .FALSE.
     811    ELSE
     812       use_top_fluxes = .TRUE.    ! because this is currently the only choice
     813    ENDIF
    806814
    807815!
    … …  
    823831          PRINT*, '    allowed with pt_surface_initial_change (/=0) = ', &
    824832                  pt_surface_initial_change
     833       ENDIF
     834       CALL local_stop
     835    ENDIF
     836
     837!
     838!-- A given temperature at the top implies Dirichlet boundary condition for
     839!-- temperature. In this case specification of a constant heat flux is
     840!-- forbidden.
     841    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
     842         top_heatflux /= 0.0 )  THEN
     843       IF ( myid == 0 )  THEN
     844          PRINT*, '+++ check_parameters:'
     845          PRINT*, '    boundary_condition: bc_pt_t = ', bc_pt_t
     846          PRINT*, '    is not allowed with constant_top_heatflux = .TRUE.'
    825847       ENDIF
    826848       CALL local_stop

  • palm/trunk/SOURCE/diffusion_e.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt
    77!
    88! Former revisions:
    … …  
    5656       REAL, DIMENSION(:,:), POINTER   ::  rif
    5757       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
    58        REAL, DIMENSION(nzb+1:nzt-1,nys:nyn) ::  dissipation, l, ll
     58       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
    5959 
    6060
    … …  
    7272             ENDIF
    7373
    74              DO  k = nzb_s_inner(j,i)+1, nzt-1
     74             DO  k = nzb_s_inner(j,i)+1, nzt
    7575!
    7676!--             Calculate the mixing length (for dissipation)
    … …  
    102102!--       Calculate the tendency terms
    103103          DO  j = nys, nyn
    104              DO  k = nzb_s_inner(j,i)+1, nzt-1
     104             DO  k = nzb_s_inner(j,i)+1, nzt
    105105
    106106                 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
    … …  
    130130          IF ( use_sgs_for_particles )  THEN
    131131             DO  j = nys, nyn
    132                 DO  k = nzb_s_inner(j,i)+1, nzt-1
     132                DO  k = nzb_s_inner(j,i)+1, nzt
    133133                   diss(k,j,i) = dissipation(k,j)
    134134                ENDDO
    … …  
    171171       REAL, DIMENSION(:,:), POINTER   ::  rif
    172172       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
    173        REAL, DIMENSION(nzb+1:nzt-1)    ::  dissipation, l, ll
     173       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
    174174
    175175
    … …  
    187187!
    188188!--    Calculate the mixing length (for dissipation)
    189        DO  k = nzb_s_inner(j,i)+1, nzt-1
     189       DO  k = nzb_s_inner(j,i)+1, nzt
    190190          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    191191          IF ( dpt_dz > 0.0 ) THEN
    … …  
    234234!--    Store dissipation if needed for calculating the sgs particle velocities
    235235       IF ( use_sgs_for_particles )  THEN
    236           DO  k = nzb_s_inner(j,i)+1, nzt-1
     236          DO  k = nzb_s_inner(j,i)+1, nzt
    237237             diss(k,j,i) = dissipation(k)
    238238          ENDDO

  • palm/trunk/SOURCE/diffusion_s.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt, fluxes can be given at top,
     7! +s_flux_t in parameter list, s_flux renamed s_flux_b
    78!
    89! Former revisions:
    … …  
    3940! Call for all grid points
    4041!------------------------------------------------------------------------------!
    41     SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux, tend )
     42    SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux_b, s_flux_t, tend )
    4243
    4344       USE control_parameters
    … …  
    5152       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt)
    5253       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    53        REAL, DIMENSION(:,:),   POINTER ::  s_flux
     54       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
    5455       REAL, DIMENSION(:,:,:), POINTER ::  kh, s
    5556
    … …  
    5859!
    5960!--          Compute horizontal diffusion
    60              DO  k = nzb_s_outer(j,i)+1, nzt-1
     61             DO  k = nzb_s_outer(j,i)+1, nzt
    6162
    6263                tend(k,j,i) = tend(k,j,i)                                     &
    … …  
    9798!
    9899!--          Compute vertical diffusion. In case that surface fluxes have been
    99 !--          presribed or computed, index k starts at nzb+2.
    100              DO  k = nzb_diff_s_inner(j,i), nzt-1
     100!--          prescribed or computed at bottom and/or top, index k starts/ends at
     101!--          nzb+2 or nzt-1, respectively.
     102             DO  k = nzb_diff_s_inner(j,i), nzt_diff
    101103
    102104                tend(k,j,i) = tend(k,j,i)                                     &
    … …  
    108110
    109111!
    110 !--          Vertical diffusion at the first computational gridpoint in &
     112!--          Vertical diffusion at the first computational gridpoint along
    111113!--          z-direction
    112114             IF ( use_surface_fluxes )  THEN
    … …  
    118120                                               * ( s(k+1,j,i)-s(k,j,i) )      &
    119121                                               * ddzu(k+1)                    &
    120                                            + s_flux(j,i)                      &
     122                                           + s_flux_b(j,i)                    &
     123                                         ) * ddzw(k)
     124
     125             ENDIF
     126
     127!
     128!--          Vertical diffusion at the last computational gridpoint along
     129!--          z-direction
     130             IF ( use_top_fluxes )  THEN
     131
     132                k = nzt
     133
     134                tend(k,j,i) = tend(k,j,i)                                     &
     135                                       + ( - s_flux_t(j,i)                    &
     136                                         - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )    &
     137                                               * ( s(k,j,i)-s(k-1,j,i) )      &
     138                                               * ddzu(k)                      &
    121139                                         ) * ddzw(k)
    122140
    … …  
    132150! Call for grid point i,j
    133151!------------------------------------------------------------------------------!
    134     SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux, tend )
     152    SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux_b, s_flux_t, &
     153                               tend )
    135154
    136155       USE control_parameters
    … …  
    144163       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt)
    145164       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    146        REAL, DIMENSION(:,:),   POINTER ::  s_flux
     165       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
    147166       REAL, DIMENSION(:,:,:), POINTER ::  kh, s
    148167
    149168!
    150169!--    Compute horizontal diffusion
    151        DO  k = nzb_s_outer(j,i)+1, nzt-1
     170       DO  k = nzb_s_outer(j,i)+1, nzt
    152171
    153172          tend(k,j,i) = tend(k,j,i)                                           &
    … …  
    188207!
    189208!--    Compute vertical diffusion. In case that surface fluxes have been
    190 !--    presribed or computed, index k starts at nzb+2.
    191        DO  k = nzb_diff_s_inner(j,i), nzt-1
     209!--    prescribed or computed at bottom and/or top, index k starts/ends at
     210!--    nzb+2 or nzt-1, respectively.
     211       DO  k = nzb_diff_s_inner(j,i), nzt_diff
    192212
    193213          tend(k,j,i) = tend(k,j,i)                                           &
    … …  
    199219
    200220!
    201 !--    Vertical diffusion at the first computational gridpoint in z-direction
     221!--    Vertical diffusion at the first computational gridpoint along z-direction
    202222       IF ( use_surface_fluxes )  THEN
    203223
    204224          k = nzb_s_inner(j,i)+1
    205225
    206           tend(k,j,i) = tend(k,j,i)                                           &
    207                                        + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )    &
    208                                                * ( s(k+1,j,i)-s(k,j,i) )      &
    209                                                * ddzu(k+1)                    &
    210                                            + s_flux(j,i)                      &
    211                                          ) * ddzw(k)
     226          tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
     227                                            * ( s(k+1,j,i)-s(k,j,i) )    &
     228                                            * ddzu(k+1)                  &
     229                                        + s_flux_b(j,i)                  &
     230                                      ) * ddzw(k)
    212231
    213232       ENDIF
    214233
     234!
     235!--    Vertical diffusion at the last computational gridpoint along z-direction
     236       IF ( use_top_fluxes )  THEN
     237
     238          k = nzt
     239
     240          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                  &
     241                                      - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )  &
     242                                            * ( s(k,j,i)-s(k-1,j,i) )    &
     243                                            * ddzu(k)                    &
     244                                      ) * ddzw(k)
     245
     246       ENDIF
     247
    215248    END SUBROUTINE diffusion_s_ij
    216249

  • palm/trunk/SOURCE/flow_statistics.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! fluxes at top modified (tswst, qswst)
    77!
    88! Former revisions:
    … …  
    347347                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    348348                                                               ) * rmask(j,i,sr)
     349             ENDDO
     350
     351             DO  k = nzb_diff_s_outer(j,i)-1, nzt_diff
    349352!
    350353!--             Heat flux w"pt"
    … …  
    417420                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
    418421                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     422                ENDIF
     423             ENDIF
     424
     425!
     426!--          Subgridscale fluxes at the top surface
     427             IF ( use_top_fluxes )  THEN
     428                sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
     429                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
     430                sums_l(nzt,58,tn) = sums_l(nzt,58,tn) + &
     431                                    0.0 * rmask(j,i,sr)           ! u"pt"
     432                sums_l(nzt,61,tn) = sums_l(nzt,61,tn) + &
     433                                    0.0 * rmask(j,i,sr)           ! v"pt"
     434                IF ( moisture )  THEN
     435                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
     436                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     437                   IF ( cloud_physics )  THEN
     438                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
     439                                          ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
     440                                          tswst(j,i) + 0.61 * pt(nzt,j,i) * &
     441                                                     qsws(j,i)            &
     442                                                              )
     443!
     444!--                   Formula does not work if ql(nzb) /= 0.0
     445                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
     446                                          qswst(j,i) * rmask(j,i,sr)
     447                   ENDIF
     448                ENDIF
     449                IF ( passive_scalar )  THEN
     450                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
     451                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    419452                ENDIF
    420453             ENDIF

  • palm/trunk/SOURCE/header.f90

    r4 r19  
    345345    ENDIF
    346346    IF ( ibc_pt_t == 0 )  THEN
    347        roben  = TRIM( roben  ) // ' pt(nzt) = pt_top'
    348     ELSE
    349        roben  = TRIM( roben  ) // ' pt(nzt) = pt(nzt-1) + dpt/dz'
     347       roben  = TRIM( roben  ) // ' pt(nzt+1) = pt_top'
     348    ELSEIF( ibc_pt_t == 1 )  THEN
     349       roben  = TRIM( roben  ) // ' pt(nzt+1) = pt(nzt)'
     350    ELSEIF( ibc_pt_t == 2 )  THEN
     351       roben  = TRIM( roben  ) // ' pt(nzt+1) = pt(nzt) + dpt/dz_ini'
    350352    ENDIF
    351353
    … …  
    404406       IF ( passive_scalar  .AND.  constant_waterflux )  THEN
    405407          WRITE ( io, 313 ) surface_waterflux
     408       ENDIF
     409    ENDIF
     410
     411    IF ( use_top_fluxes )  THEN
     412       WRITE ( io, 304 )
     413       IF ( constant_top_heatflux )  THEN
     414          WRITE ( io, 306 )  top_heatflux
     415       ENDIF
     416       IF ( moisture  .OR.  passive_scalar )  THEN
     417          WRITE ( io, 315 )
    406418       ENDIF
    407419    ENDIF
    … …  
    12201232             ' B. bound.: ',A/ &
    12211233             ' T. bound.: ',A)
    1222 303 FORMAT (/' Surface fluxes are used in diffusion terms at k=1')
    1223 305 FORMAT (//'    Prandtl-Layer between surface and first computational ', &
    1224                'u,v-level:'// &
     1234303 FORMAT (/' Bottom surface fluxes are used in diffusion terms at k=1')
     1235304 FORMAT (/' Top surface fluxes are used in diffusion terms at k=nzt')
     1236305 FORMAT (//'    Prandtl-Layer between bottom surface and first ', &
     1237               'computational u,v-level:'// &
    12251238             '       zp = ',F6.2,' m   z0 = ',F6.4,' m   kappa = ',F4.2/ &
    12261239             '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
    … …  
    12341247313 FORMAT ('       Predefined constant scalar flux: ',E10.3,' kg/(m**2 s)')
    12351248314 FORMAT ('       Predefined scalar value at the surface')
     1249315 FORMAT ('       Humidity / scalar flux at top surface is 0.0')
    12361250317 FORMAT (//' Lateral boundaries:'/ &
    12371251            '       left/right:  ',A/    &

  • palm/trunk/SOURCE/impact_of_latent_heat.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt
    77!
    88! Former revisions:
    … …  
    5353       DO  i = nxl, nxr
    5454          DO  j = nys, nyn
    55              DO  k = nzb_2d(j,i)+1, nzt-1   
     55             DO  k = nzb_2d(j,i)+1, nzt
    5656
    5757                IF ( ql(k,j,i) > ql_crit )  THEN
    … …  
    8787
    8888
    89        DO  k = nzb_2d(j,i)+1, nzt-1   
     89       DO  k = nzb_2d(j,i)+1, nzt
    9090
    9191          IF ( ql(k,j,i) > ql_crit )  THEN

  • palm/trunk/SOURCE/init_3d_model.f90

    r4 r19  
    77! Actual revisions:
    88! -----------------
    9 !
     9! +handling of top fluxes
    1010!
    1111! Former revisions:
    … …  
    7777    ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
    7878
    79     ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1),  &
    80               ts(nys-1:nyn+1,nxl-1:nxr+1), us(nys-1:nyn+1,nxl-1:nxr+1),        &
    81               usws_1(nys-1:nyn+1,nxl-1:nxr+1), vsws_1(nys-1:nyn+1,nxl-1:nxr+1),&
    82               z0(nys-1:nyn+1,nxl-1:nxr+1) )
     79    ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
     80              ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
     81              us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
     82              vsws_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
    8383
    8484    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    8585!
    8686!--    Leapfrog scheme needs two timelevels of diffusion quantities
    87        ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),  &
    88                  shf_2(nys-1:nyn+1,nxl-1:nxr+1),  &
    89                  usws_2(nys-1:nyn+1,nxl-1:nxr+1), &
     87       ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),   &
     88                 shf_2(nys-1:nyn+1,nxl-1:nxr+1),   &
     89                 tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
     90                 usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
    9091                 vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
    9192    ENDIF
    … …  
    121122!--    2D-moisture/scalar arrays
    122123       ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1),     &
    123                   qsws_1(nys-1:nyn+1,nxl-1:nxr+1) )
     124                  qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
     125                  qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
    124126
    125127       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    126           ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
     128          ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
     129                    qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
    127130       ENDIF
    128131!
    … …  
    177180    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    178181
    179        rif_m  => rif_1;   rif  => rif_2
    180        shf_m  => shf_1;   shf  => shf_2
    181        usws_m => usws_1;  usws => usws_2
    182        vsws_m => vsws_1;  vsws => vsws_2
     182       rif_m   => rif_1;    rif   => rif_2
     183       shf_m   => shf_1;    shf   => shf_2
     184       tswst_m => tswst_1;  tswst => tswst_2
     185       usws_m  => usws_1;   usws  => usws_2
     186       vsws_m  => vsws_1;   vsws  => vsws_2
    183187       e_m  => e_1;   e  => e_2;   e_p  => e_3;   te_m  => e_3
    184188       kh_m => kh_1;  kh => kh_2
    … …  
    190194
    191195       IF ( moisture  .OR.  passive_scalar )  THEN
    192           qsws_m => qsws_1;  qsws => qsws_2
     196          qsws_m  => qsws_1;   qsws  => qsws_2
     197          qswst_m => qswst_1;  qswst => qswst_2
    193198          q_m    => q_1;     q    => q_2;     q_p => q_3;     tq_m => q_3
    194199          IF ( moisture )        vpt_m  => vpt_1;   vpt  => vpt_2
    … …  
    202207    ELSE
    203208
    204        rif  => rif_1
    205        shf  => shf_1
    206        usws => usws_1
    207        vsws => vsws_1
    208        e    => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
    209        kh   => kh_1
    210        km   => km_1
    211        pt   => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
    212        u    => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
    213        v    => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
    214        w    => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
     209       rif   => rif_1
     210       shf   => shf_1
     211       tswst => tswst_1
     212       usws  => usws_1
     213       vsws  => vsws_1
     214       e     => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
     215       kh    => kh_1
     216       km    => km_1
     217       pt    => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
     218       u     => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
     219       v     => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
     220       w     => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
    215221
    216222       IF ( moisture  .OR.  passive_scalar )  THEN
    217223          qsws   => qsws_1
     224          qswst  => qswst_1
    218225          q      => q_1;     q_p  => q_2;     tq_m => q_3;    q_m => q_3
    219226          IF ( moisture )        vpt  => vpt_1
    … …  
    453460
    454461!
    455 !--    Initialize surface fluxes
     462!--    Initialize fluxes at bottom surface
    456463       IF ( use_surface_fluxes )  THEN
    457464
    … …  
    486493             ENDIF
    487494          ENDIF
     495
     496       ENDIF
     497
     498!
     499!--    Initialize fluxes at top surface
     500!--    Currently, only the heatflux can be prescribed. The latent flux is
     501!--    zeri in this case!
     502       IF ( use_top_fluxes )  THEN
     503
     504          IF ( constant_top_heatflux )  THEN
     505!
     506!--          Heat flux is prescribed
     507             tswst = top_heatflux
     508             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
     509
     510             IF ( moisture  .OR.  passive_scalar )  THEN
     511                qswst = 0.0
     512                IF ( ASSOCIATED( qswst_m ) )  qswst_m = qswst
     513             ENDIF
     514         ENDIF
    488515
    489516       ENDIF

  • palm/trunk/SOURCE/init_grid.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Setting of nzt_diff
    77!
    88! Former revisions:
    … …  
    202202
    203203!
    204 !-- Define vertical gridpoint from which on the usual finite difference
     204!-- Define vertical gridpoint from (or to) which on the usual finite difference
    205205!-- form (which does not use surface fluxes) is applied
    206206    IF ( prandtl_layer  .OR.  use_surface_fluxes )  THEN
    … …  
    208208    ELSE
    209209       nzb_diff = nzb + 1
     210    ENDIF
     211    IF ( use_top_fluxes )  THEN
     212       nzt_diff = nzt - 1
     213    ELSE
     214       nzt_diff = nzt
    210215    ENDIF
    211216

  • palm/trunk/SOURCE/init_pt_anomaly.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt
    77!
    88! Former revisions:
    … …  
    4747    DO  i = nxl, nxr
    4848       DO  j = nys, nyn
    49           DO  k = nzb+1, nzt-1
     49          DO  k = nzb+1, nzt
    5050             x = ( i - ic ) * dx
    5151             y = ( j - jc ) * dy

  • palm/trunk/SOURCE/modules.f90

    r4 r19  
    55! Actual revisions:
    66! -----------------
    7 !
     7! +constant_top_heatflux, top_heatflux, use_top_fluxes, +arrays for top fluxes,
     8! +nzt_diff, default of bc_pt_t renamed "initial_gradient"
     9! Bugfix: p is not a pointer
    810!
    911! Former revisions:
    … …  
    7577
    7678    REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
    77           qsws_1, qsws_2, rif_1, rif_2, shf_1, shf_2, usws_1, usws_2,          &
    78           vsws_1, vsws_2
     79          qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, shf_1, shf_2,        &
     80          tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2
    7981
    8082    REAL, DIMENSION(:,:), POINTER ::                                           &
    81           qsws, qsws_m, rif, rif_m, shf, shf_m, usws, usws_m, vsws, vsws_m
     83          qsws, qsws_m, qswst, qswst_m, rif, rif_m, shf, shf_m, tswst,         &
     84          tswst_m, usws, usws_m, vsws, vsws_m
    8285
    8386    REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
    … …  
    8891
    8992    REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
    90           e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, pt_1, pt_2, pt_3, q_1, q_2,   &
    91           q_3, ql_1, ql_2, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, vpt_2, w_1,    &
    92           w_2, w_3
     93          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, pt_1, pt_2, pt_3, q_1,     &
     94          q_2, q_3, ql_1, ql_2, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, vpt_2,    &
     95          w_1, w_2, w_3
    9396
    9497    REAL, DIMENSION(:,:,:), POINTER ::                                         &
    95           e, e_m, e_p, kh, kh_m, km, km_m, p, pt, pt_m, pt_p, q, q_m, q_p, ql, &
     98          e, e_m, e_p, kh, kh_m, km, km_m, pt, pt_m, pt_p, q, q_m, q_p, ql,    &
    9699          ql_c, te_m, tpt_m, tq_m, tu_m, tv_m, tw_m, u, u_m, u_p, v, v_m, v_p, &
    97100          vpt, vpt_m, w, w_m, w_p
    … …  
    201204                            scalar_advec = 'pw-scheme'
    202205    CHARACTER (LEN=20)  ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
    203                             bc_ns = 'cyclic', &
    204                             bc_p_b = 'neumann', bc_p_t = 'dirichlet', &
    205                             bc_pt_b = 'dirichlet', bc_pt_t = 'neumann', &
     206                            bc_ns = 'cyclic', bc_p_b = 'neumann', &
     207                            bc_p_t = 'dirichlet', bc_pt_b = 'dirichlet', &
     208                            bc_pt_t = 'initial_gradient', &
    206209                            bc_q_b = 'dirichlet', bc_q_t = 'neumann', &
    207210                            bc_s_b = 'dirichlet', bc_s_t = 'neumann', &
    … …  
    263266                cloud_droplets = .FALSE., cloud_physics = .FALSE., &
    264267                conserve_volume_flow = .FALSE., constant_diffusion = .FALSE., &
    265                 constant_heatflux = .TRUE., constant_waterflux = .TRUE., &
    266                 create_disturbances = .TRUE., cut_spline_overshoot = .TRUE., &
     268                constant_heatflux = .TRUE., constant_top_heatflux = .TRUE., &
     269                constant_waterflux = .TRUE., create_disturbances = .TRUE., &
     270                cut_spline_overshoot = .TRUE., &
    267271                data_output_2d_on_each_pe = .TRUE., do2d_at_begin = .FALSE., &
    268272                do3d_at_begin = .FALSE., do3d_compress = .FALSE., &
    … …  
    270274                disturbance_created = .FALSE., &
    271275                first_call_advec_particles = .TRUE., &
    272                 force_print_header = .FALSE., galilei_transformation = .FALSE., &
     276                force_print_header = .FALSE., galilei_transformation = .FALSE.,&
    273277                inflow_l = .FALSE., inflow_n = .FALSE., inflow_r = .FALSE., &
    274278                inflow_s = .FALSE., iso2d_output = .FALSE., &
    … …  
    281285                random_heatflux = .FALSE., run_control_header = .FALSE., &
    282286                sloping_surface = .FALSE., stop_dt = .FALSE., &
    283                 terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE., &
    284                 use_surface_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., &
     287                terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,&
     288                use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
     289                use_ug_for_galilei_tr = .TRUE., &
    285290                use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE.
    286291
    … …  
    336341             time_do_sla = 0.0, time_dvrp = 0.0, time_prel = 0.0, &
    337342             time_restart = 9999999.9, time_run_control = 0.0, &
    338              ug_surface = 0.0, u_gtrans = 0.0, ups_limit_e = 0.0, &
    339              ups_limit_pt = 0.0, ups_limit_u = 0.0, ups_limit_v = 0.0, &
    340              ups_limit_w = 0.0, vg_surface = 0.0, v_gtrans = 0.0, &
    341              wall_adjustment_factor = 1.8, z_max_do1d = -1.0, &
     343             top_heatflux = 9999999.9, ug_surface = 0.0, u_gtrans = 0.0, &
     344             ups_limit_e = 0.0, ups_limit_pt = 0.0, ups_limit_u = 0.0, &
     345             ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, &
     346             v_gtrans = 0.0, wall_adjustment_factor = 1.8, z_max_do1d = -1.0, &
    342347             z_max_do1d_normalized = -1.0, z_max_do2d = -1.0
    343348
    … …  
    500505    INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxr, nxra, nny, ny = 0, nya,  &
    501506                nyn, nyna, nys, nnz, nz = 0, nza, nzb, nzb_diff, nzt, nzta,    &
    502                 uxrp = 0, vynp = 0
     507                nzt_diff, uxrp = 0, vynp = 0
    503508
    504509    INTEGER, DIMENSION(:), ALLOCATABLE ::                                      &

  • palm/trunk/SOURCE/parin.f90

    r7 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +top_heatflux in inipar
    77!
    88! Former revisions:
    … …  
    7171                       s_surface, s_surface_initial_change, &
    7272                       s_vertical_gradient, s_vertical_gradient_level, &
    73                        timestep_scheme, topography, ug_surface, &
     73                       top_heatflux, timestep_scheme, topography, ug_surface, &
    7474                       ug_vertical_gradient, ug_vertical_gradient_level, &
    7575                       ups_limit_e, ups_limit_pt, ups_limit_u, ups_limit_v, &

  • palm/trunk/SOURCE/production_e.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation extended for gridpoint nzt, extended for given temperature /
     7! humidity fluxes at the top
    78!
    89! Former revisions:
    … …  
    6768
    6869          DO  j = nys, nyn
    69              DO  k = nzb_diff_s_outer(j,i), nzt-1
     70             DO  k = nzb_diff_s_outer(j,i), nzt
    7071
    7172                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    … …  
    327328             DO  j = nys, nyn
    328329
    329                 DO  k = nzb_diff_s_inner(j,i), nzt-1
     330                DO  k = nzb_diff_s_inner(j,i), nzt_diff
    330331                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
    331332                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
    332333                ENDDO
     334
    333335                IF ( use_surface_fluxes )  THEN
    334336                   k = nzb_diff_s_inner(j,i)-1
    … …  
    336338                ENDIF
    337339
     340                IF ( use_top_fluxes )  THEN
     341                   k = nzt
     342                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
     343                ENDIF
     344
    338345             ENDDO
    339346
    … …  
    342349             DO  j = nys, nyn
    343350
    344                 DO  k = nzb_diff_s_inner(j,i), nzt-1
     351                DO  k = nzb_diff_s_inner(j,i), nzt_diff
    345352
    346353                   IF ( .NOT. cloud_physics )  THEN
    … …  
    375382                DO  j = nys, nyn
    376383
    377                    k = nzb_diff_s_inner(j,i)-1
     384                   k = nzb_diff_s_inner(j,i)
    378385
    379386                   IF ( .NOT. cloud_physics )  THEN
    … …  
    402409             ENDIF
    403410
     411             IF ( use_top_fluxes )  THEN
     412
     413                DO  j = nys, nyn
     414
     415                   k = nzt
     416
     417                   IF ( .NOT. cloud_physics )  THEN
     418                      k1 = 1.0 + 0.61 * q(k,j,i)
     419                      k2 = 0.61 * pt(k,j,i)
     420                   ELSE
     421                      IF ( ql(k,j,i) == 0.0 )  THEN
     422                         k1 = 1.0 + 0.61 * q(k,j,i)
     423                         k2 = 0.61 * pt(k,j,i)
     424                      ELSE
     425                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
     426                         temp  = theta * t_d_pt(k)
     427                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
     428                                    ( q(k,j,i) - ql(k,j,i) ) *          &
     429                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
     430                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     431                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
     432                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     433                      ENDIF
     434                   ENDIF
     435
     436                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
     437                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
     438                ENDDO
     439
     440             ENDIF
     441
    404442          ENDIF
    405443
    … …  
    430468!
    431469!--    Calculate TKE production by shear
    432        DO  k = nzb_diff_s_outer(j,i), nzt-1
     470       DO  k = nzb_diff_s_outer(j,i), nzt
    433471
    434472          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    … …  
    657695       IF ( .NOT. moisture )  THEN
    658696
    659           DO  k = nzb_diff_s_inner(j,i), nzt-1
     697          DO  k = nzb_diff_s_inner(j,i), nzt_diff
    660698             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
    661699                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
    662700          ENDDO
     701
    663702          IF ( use_surface_fluxes )  THEN
    664703             k = nzb_diff_s_inner(j,i)-1
    … …  
    666705          ENDIF
    667706
     707          IF ( use_top_fluxes )  THEN
     708             k = nzt
     709             tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
     710          ENDIF
     711
    668712       ELSE
    669713
    670           DO  k = nzb_diff_s_inner(j,i), nzt-1
     714          DO  k = nzb_diff_s_inner(j,i), nzt_diff
    671715
    672716             IF ( .NOT. cloud_physics )  THEN
    … …  
    694738                                      ) * dd2zu(k)
    695739          ENDDO
     740
    696741          IF ( use_surface_fluxes )  THEN
    697742             k = nzb_diff_s_inner(j,i)-1
    … …  
    720765          ENDIF
    721766
     767          IF ( use_top_fluxes )  THEN
     768             k = nzt
     769
     770             IF ( .NOT. cloud_physics ) THEN
     771                k1 = 1.0 + 0.61 * q(k,j,i)
     772                k2 = 0.61 * pt(k,j,i)
     773             ELSE
     774                IF ( ql(k,j,i) == 0.0 )  THEN
     775                   k1 = 1.0 + 0.61 * q(k,j,i)
     776                   k2 = 0.61 * pt(k,j,i)
     777                ELSE
     778                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
     779                   temp  = theta * t_d_pt(k)
     780                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
     781                              ( q(k,j,i) - ql(k,j,i) ) *          &
     782                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
     783                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     784                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
     785                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     786                ENDIF
     787             ENDIF
     788
     789             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
     790                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
     791          ENDIF
     792
    722793       ENDIF
    723794

  • palm/trunk/SOURCE/prognostic_equations.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation of e, q, and pt extended for gridpoint nzt,
     7! handling of given temperature/humidity/scalar fluxes at top surface
    78!
    89! Former revisions:
    … …  
    2223! Description:
    2324! ------------
    24 ! Solving the prognostic equations and advecting particles.
     25! Solving the prognostic equations.
    2526!------------------------------------------------------------------------------!
    2627
    … …  
    333334!--       Tendency terms
    334335          IF ( scalar_advec == 'bc-scheme' )  THEN
    335              CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tend )
     336             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
    336337          ELSE
    337338             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
    … …  
    346347             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
    347348             THEN
    348                 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, tend )
     349                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
     350                                  tswst_m, tend )
    349351             ELSE
    350                 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tend )
     352                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
    351353             ENDIF
    352354          ENDIF
    … …  
    368370!
    369371!--       Prognostic equation for potential temperature
    370           DO  k = nzb_s_inner(j,i)+1, nzt-1
     372          DO  k = nzb_s_inner(j,i)+1, nzt
    371373             pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
    372374                           dt_3d * (                                           &
    … …  
    380382          IF ( timestep_scheme(1:5) == 'runge' )  THEN
    381383             IF ( intermediate_timestep_count == 1 )  THEN
    382                 DO  k = nzb_s_inner(j,i)+1, nzt-1
     384                DO  k = nzb_s_inner(j,i)+1, nzt
    383385                   tpt_m(k,j,i) = tend(k,j,i)
    384386                ENDDO
    385387             ELSEIF ( intermediate_timestep_count < &
    386388                      intermediate_timestep_count_max )  THEN
    387                 DO  k = nzb_s_inner(j,i)+1, nzt-1
     389                DO  k = nzb_s_inner(j,i)+1, nzt
    388390                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
    389391                ENDDO
    … …  
    429431!--          Tendency-terms
    430432             IF ( scalar_advec == 'bc-scheme' )  THEN
    431                 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )
     433                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, tend )
    432434             ELSE
    433435                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) THEN
    … …  
    443445                THEN
    444446                   CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
     447                                     qswst_m, tend )
     448                ELSE
     449                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
    445450                                     tend )
    446                 ELSE
    447                    CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )
    448451                ENDIF
    449452             ENDIF
    … …  
    459462!
    460463!--          Prognostic equation for total water content / scalar
    461              DO  k = nzb_s_inner(j,i)+1, nzt-1
     464             DO  k = nzb_s_inner(j,i)+1, nzt
    462465                q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) +       &
    463466                             dt_3d * (                                         &
    … …  
    471474             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    472475                IF ( intermediate_timestep_count == 1 )  THEN
    473                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     476                   DO  k = nzb_s_inner(j,i)+1, nzt
    474477                      tq_m(k,j,i) = tend(k,j,i)
    475478                   ENDDO
    476479                ELSEIF ( intermediate_timestep_count < &
    477480                         intermediate_timestep_count_max )  THEN
    478                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     481                   DO  k = nzb_s_inner(j,i)+1, nzt
    479482                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
    480483                   ENDDO
    … …  
    577580!--          reasons in the course of the integration. In such cases the old TKE
    578581!--          value is reduced by 90%.
    579              DO  k = nzb_s_inner(j,i)+1, nzt-1
     582             DO  k = nzb_s_inner(j,i)+1, nzt
    580583                e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) +       &
    581584                             dt_3d * (                                         &
    … …  
    589592             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    590593                IF ( intermediate_timestep_count == 1 )  THEN
    591                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     594                   DO  k = nzb_s_inner(j,i)+1, nzt
    592595                      te_m(k,j,i) = tend(k,j,i)
    593596                   ENDDO
    594597                ELSEIF ( intermediate_timestep_count < &
    595598                         intermediate_timestep_count_max )  THEN
    596                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     599                   DO  k = nzb_s_inner(j,i)+1, nzt
    597600                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
    598601                   ENDDO
    … …  
    809812             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
    810813             THEN
    811                 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, tend )
     814                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
     815                                  tswst_m, tend )
    812816             ELSE
    813                 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tend )
     817                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
    814818             ENDIF
    815819
    … …  
    830834!
    831835!--          Prognostic equation for potential temperature
    832              DO  k = nzb_s_inner(j,i)+1, nzt-1
     836             DO  k = nzb_s_inner(j,i)+1, nzt
    833837                pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
    834838                              dt_3d * (                                        &
    … …  
    842846             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    843847                IF ( intermediate_timestep_count == 1 )  THEN
    844                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     848                   DO  k = nzb_s_inner(j,i)+1, nzt
    845849                      tpt_m(k,j,i) = tend(k,j,i)
    846850                   ENDDO
    847851                ELSEIF ( intermediate_timestep_count < &
    848852                         intermediate_timestep_count_max )  THEN
    849                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     853                   DO  k = nzb_s_inner(j,i)+1, nzt
    850854                      tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
    851855                                      5.3125 * tpt_m(k,j,i)
    … …  
    870874                IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
    871875                THEN
    872                    CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, tend )
     876                   CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
     877                                     qswst_m, tend )
    873878                ELSE
    874                    CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )
     879                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
     880                                     tend )
    875881                ENDIF
    876882       
    … …  
    885891!
    886892!--             Prognostic equation for total water content / scalar
    887                 DO  k = nzb_s_inner(j,i)+1, nzt-1
     893                DO  k = nzb_s_inner(j,i)+1, nzt
    888894                   q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
    889895                                dt_3d * (                                      &
    … …  
    897903                IF ( timestep_scheme(1:5) == 'runge' )  THEN
    898904                   IF ( intermediate_timestep_count == 1 )  THEN
    899                       DO  k = nzb_s_inner(j,i)+1, nzt-1
     905                      DO  k = nzb_s_inner(j,i)+1, nzt
    900906                         tq_m(k,j,i) = tend(k,j,i)
    901907                      ENDDO
    902908                   ELSEIF ( intermediate_timestep_count < &
    903909                            intermediate_timestep_count_max )  THEN
    904                       DO  k = nzb_s_inner(j,i)+1, nzt-1
     910                      DO  k = nzb_s_inner(j,i)+1, nzt
    905911                         tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
    906912                                        5.3125 * tq_m(k,j,i)
    … …  
    951957!--             reasons in the course of the integration. In such cases the old
    952958!--             TKE value is reduced by 90%.
    953                 DO  k = nzb_s_inner(j,i)+1, nzt-1
     959                DO  k = nzb_s_inner(j,i)+1, nzt
    954960                   e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
    955961                                dt_3d * (                                      &
    … …  
    963969                IF ( timestep_scheme(1:5) == 'runge' )  THEN
    964970                   IF ( intermediate_timestep_count == 1 )  THEN
    965                       DO  k = nzb_s_inner(j,i)+1, nzt-1
     971                      DO  k = nzb_s_inner(j,i)+1, nzt
    966972                         te_m(k,j,i) = tend(k,j,i)
    967973                      ENDDO
    968974                   ELSEIF ( intermediate_timestep_count < &
    969975                            intermediate_timestep_count_max )  THEN
    970                       DO  k = nzb_s_inner(j,i)+1, nzt-1
     976                      DO  k = nzb_s_inner(j,i)+1, nzt
    971977                         te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
    972978                                        5.3125 * te_m(k,j,i)
    … …  
    12521258!-- pt-tendency terms with no communication
    12531259    IF ( scalar_advec == 'bc-scheme' )  THEN
    1254        CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tend )
     1260       CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
    12551261    ELSE
    12561262       IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
    … …  
    12641270       ENDIF
    12651271       IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    1266           CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tend )
     1272          CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, tend )
    12671273       ELSE
    1268           CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tend )
     1274          CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
    12691275       ENDIF
    12701276    ENDIF
    … …  
    12881294    DO  i = nxl, nxr
    12891295       DO  j = nys, nyn
    1290           DO  k = nzb_s_inner(j,i)+1, nzt-1
     1296          DO  k = nzb_s_inner(j,i)+1, nzt
    12911297             pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +       &
    12921298                           dt_3d * (                                           &
    … …  
    13041310          DO  i = nxl, nxr
    13051311             DO  j = nys, nyn
    1306                 DO  k = nzb_s_inner(j,i)+1, nzt-1
     1312                DO  k = nzb_s_inner(j,i)+1, nzt
    13071313                   tpt_m(k,j,i) = tend(k,j,i)
    13081314                ENDDO
    … …  
    13131319          DO  i = nxl, nxr
    13141320             DO  j = nys, nyn
    1315                 DO  k = nzb_s_inner(j,i)+1, nzt-1
     1321                DO  k = nzb_s_inner(j,i)+1, nzt
    13161322                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
    13171323                ENDDO
    … …  
    13521358!--    Scalar/q-tendency terms with no communication
    13531359       IF ( scalar_advec == 'bc-scheme' )  THEN
    1354           CALL diffusion_s( ddzu, ddzw, kh, q, qsws, tend )
     1360          CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
    13551361       ELSE
    13561362          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
    … …  
    13641370          ENDIF
    13651371          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    1366              CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, tend )
     1372             CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, tend )
    13671373          ELSE
    1368              CALL diffusion_s( ddzu, ddzw, kh, q, qsws, tend )
     1374             CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
    13691375          ENDIF
    13701376       ENDIF
    … …  
    13821388       DO  i = nxl, nxr
    13831389          DO  j = nys, nyn
    1384              DO  k = nzb_s_inner(j,i)+1, nzt-1
     1390             DO  k = nzb_s_inner(j,i)+1, nzt
    13851391                q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) +       &
    13861392                             dt_3d * (                                         &
    … …  
    13981404             DO  i = nxl, nxr
    13991405                DO  j = nys, nyn
    1400                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     1406                   DO  k = nzb_s_inner(j,i)+1, nzt
    14011407                      tq_m(k,j,i) = tend(k,j,i)
    14021408                   ENDDO
    … …  
    14071413             DO  i = nxl, nxr
    14081414                DO  j = nys, nyn
    1409                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     1415                   DO  k = nzb_s_inner(j,i)+1, nzt
    14101416                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
    14111417                   ENDDO
    … …  
    15031509       DO  i = nxl, nxr
    15041510          DO  j = nys, nyn
    1505              DO  k = nzb_s_inner(j,i)+1, nzt-1
     1511             DO  k = nzb_s_inner(j,i)+1, nzt
    15061512                e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) +       &
    15071513                             dt_3d * (                                         &
    … …  
    15191525             DO  i = nxl, nxr
    15201526                DO  j = nys, nyn
    1521                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     1527                   DO  k = nzb_s_inner(j,i)+1, nzt
    15221528                      te_m(k,j,i) = tend(k,j,i)
    15231529                   ENDDO
    … …  
    15281534             DO  i = nxl, nxr
    15291535                DO  j = nys, nyn
    1530                    DO  k = nzb_s_inner(j,i)+1, nzt-1
     1536                   DO  k = nzb_s_inner(j,i)+1, nzt
    15311537                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
    15321538                   ENDDO

  • palm/trunk/SOURCE/read_3d_binary.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +qswst, qswst_m, tswst, tswst_m
    77!
    88! Former revisions:
    … …  
    249249          CASE ( 'qsws_m' )
    250250             READ ( 13 )  qsws_m
     251          CASE ( 'qswst' )
     252             READ ( 13 )  qswst
     253          CASE ( 'qswst_m' )
     254             READ ( 13 )  qswst_m
    251255          CASE ( 'qv_av' )
    252256             ALLOCATE( qv_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    … …  
    266270          CASE ( 'shf_m' )
    267271             READ ( 13 )  shf_m
     272          CASE ( 'tswst' )
     273             READ ( 13 )  tswst
     274          CASE ( 'tswst_m' )
     275             READ ( 13 )  tswst_m
    268276          CASE ( 'spectrum_x' )
    269277             READ ( 13 )  spectrum_x

  • palm/trunk/SOURCE/read_var_list.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +top_heatflux
    77!
    88! Former revisions:
    … …  
    340340          CASE ( 'topography' )
    341341             READ ( 13 )  topography
     342          CASE ( 'top_heatflux' )
     343             READ ( 13 )  top_heatflux
    342344          CASE ( 'tsc' )
    343345             READ ( 13 )  tsc

  • palm/trunk/SOURCE/spline_z.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Boundary condition for pt at top adjusted
    77!
    88! Former revisions:
    … …  
    156156!
    157157!--       Top boundary for temperature
    158           IF ( ibc_pt_t == 1 )  THEN
    159              vad(nzt,:)   = vad(nzt-1,:) + bc_pt_t_val * dz_spline(nzt)
     158          IF ( ibc_pt_t == 0 )  THEN
     159             vad(nzt+1,:) = pt(nzt+1,nys:nyn,i)
     160          ELSEIF ( ibc_pt_t == 1 )  THEN
     161             vad(nzt+1,:) = vad(nzt,:)
     162          ELSEIF ( ibc_pt_t == 2 )  THEN
    160163             vad(nzt+1,:) = vad(nzt,:)   + bc_pt_t_val * dz_spline(nzt+1)
    161           ELSE
    162              vad(nzt,:)   = pt(nzt,nys:nyn,i)
    163              vad(nzt+1,:) = pt(nzt+1,nys:nyn,i)
    164164          ENDIF
    165165

  • palm/trunk/SOURCE/swap_timelevel.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Swaping of top fluxes
    77!
    88! Former revisions:
    … …  
    122122                   rif_m  => rif_1;   rif  => rif_2
    123123                ENDIF
     124                IF ( use_top_fluxes )  THEN
     125                   tswst_m => tswst_1;  tswst => tswst_2
     126                   IF ( moisture  .OR.  passive_scalar )  THEN
     127                      qswst_m => qswst_1;  qswst => qswst_2
     128                   ENDIF
     129                ENDIF
    124130             ENDIF
    125131
    … …  
    161167                   rif_m  => rif_2;   rif  => rif_1
    162168                ENDIF
     169                IF ( use_top_fluxes )  THEN
     170                   tswst_m  => tswst_2;  tswst => tswst_1
     171                   IF ( moisture  .OR.  passive_scalar )  THEN
     172                      qswst_m => qswst_2;  qswst => qswst_1
     173                   ENDIF
     174                ENDIF
    163175             ENDIF
    164176

  • palm/trunk/SOURCE/write_3d_binary.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +qswst, qswst_m, tswst, tswst_m
    77!
    88! Former revisions:
    … …  
    111111          WRITE ( 14 )  'qsws_m              ';  WRITE ( 14 ) qsws_m
    112112       ENDIF
     113       WRITE ( 14 )  'qswst               ';  WRITE ( 14 ) qswst
     114       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     115          WRITE ( 14 )  'qswst_m             ';  WRITE ( 14 ) qswst_m
     116       ENDIF
    113117    ENDIF
    114118    IF ( ALLOCATED( ql_c_av ) )  THEN
    … …  
    140144    IF ( ALLOCATED( ts_av ) )  THEN
    141145       WRITE ( 14 )  'ts_av               ';  WRITE ( 14 )  ts_av
     146    ENDIF
     147    WRITE ( 14 )  'tswst               ';  WRITE ( 14 )  tswst
     148    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     149       WRITE ( 14 )  'tswst_m             ';  WRITE ( 14 )  tswst_m
    142150    ENDIF
    143151    WRITE ( 14 )  'u                   ';  WRITE ( 14 )  u

  • palm/trunk/SOURCE/write_var_list.f90

    r4 r19  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +top_heatflux
    77!
    88! Former revisions:
    … …  
    306306    WRITE ( 14 )  'topography                    '
    307307    WRITE ( 14 )  topography
     308    WRITE ( 14 )  'top_heatflux                  '
     309    WRITE ( 14 )  top_heatflux
    308310    WRITE ( 14 )  'tsc                           '
    309311    WRITE ( 14 )  tsc
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |