Changeset 2508 for palm/trunk/SOURCE


Ignore:
Timestamp:
Oct 2, 2017 8:57:09 AM (4 years ago)
Author:
suehring
Message:

Bugfixes in SGS-TKE buoyancy production; revised initialization of vertical-gradient levels in case of ocean runs

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r2422 r2508  
    2525! -----------------
    2626! $Id$
     27! Bugfix, change default value of vertical_gradient level in order to consider
     28! also ocean runs
     29!
     30! 2422 2017-09-08 08:25:41Z raasch
    2731! error message in case of missing "restart" file activation string
    2832!
     
    40314035!
    40324036!--    In case of no given gradients, choose zero gradient conditions
    4033        IF ( vertical_gradient_level(1) == -1.0_wp )  THEN
     4037       IF ( vertical_gradient_level(1) == -999999.9_wp )  THEN
    40344038          vertical_gradient_level(1) = 0.0_wp
    40354039       ENDIF
  • palm/trunk/SOURCE/modules.f90

    r2499 r2508  
    2525! -----------------
    2626! $Id$
     27! Change default value for pt/q/s/sa_vertical_gradient_level
     28!
     29! 2499 2017-09-22 16:47:58Z kanani
    2730! Default changed to fft_method = 'temperton-algorithm'
    2831!
     
    14101413    REAL(wp) ::  mask_scale(3)                                     !< collective array for mask_scale_x/y/z
    14111414    REAL(wp) ::  pt_vertical_gradient(10) = 0.0_wp                 !< namelist parameter
    1412     REAL(wp) ::  pt_vertical_gradient_level(10) = -1.0_wp          !< namelist parameter
     1415    REAL(wp) ::  pt_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
    14131416    REAL(wp) ::  q_vertical_gradient(10) = 0.0_wp                  !< namelist parameter
    1414     REAL(wp) ::  q_vertical_gradient_level(10) = -1.0_wp           !< namelist parameter
     1417    REAL(wp) ::  q_vertical_gradient_level(10) = -999999.9_wp    !< namelist parameter
    14151418    REAL(wp) ::  s_vertical_gradient(10) = 0.0_wp                  !< namelist parameter
    1416     REAL(wp) ::  s_vertical_gradient_level(10) = -1.0_wp           !< namelist parameter
     1419    REAL(wp) ::  s_vertical_gradient_level(10) = -999999.9_wp    !< namelist parameter
    14171420    REAL(wp) ::  sa_vertical_gradient(10) = 0.0_wp                 !< namelist parameter
    1418     REAL(wp) ::  sa_vertical_gradient_level(10) = -1.0_wp          !< namelist parameter
     1421    REAL(wp) ::  sa_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
    14191422    REAL(wp) ::  skip_time_domask(max_masks) = 9999999.9_wp        !< namelist parameter
    14201423    REAL(wp) ::  threshold(20) = 0.0_wp                            !< namelist parameter
  • palm/trunk/SOURCE/production_e.f90

    r2478 r2508  
    2525! -----------------
    2626! $Id$
     27! - Bugfix in buoyancy production term, wrong base state was set.
     28! - Consider use_single_reference_value case if humidity is used.
     29! - In case of use_top_fluxes, use correct inverse density at model top
     30! - Consider use_surface_fluxes and use_top_fluxes in ocean case 
     31!
     32! 2478 2017-09-18 13:37:24Z suehring
    2733! Bugfix, consider case where no constant-flux layer and no surfaces fluxes
    2834! are used
     
    152158                  humidity, kappa, neutral, ocean, pt_reference,               &
    153159                  rho_reference, use_single_reference_value,                   &
    154                   use_surface_fluxes, use_top_fluxes
     160                  use_surface_fluxes, use_top_fluxes, vpt_reference
    155161
    156162       USE grid_variables,                                                     &
     
    372378                DO  m = surf_s, surf_e
    373379                   k = surf_lsm_h%k(m)
    374 !
    375 !--                Please note, actually, an interpolation of u_0 and v_0
    376 !--                onto the grid center would be required. However, this
    377 !--                would require several data transfers between 2D-grid and
    378 !--                wall type. The effect of this missing interpolation is
    379 !--                negligible. (See also production_e_init).
     380
    380381                   dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
    381382                   dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
     
    388389                DO  m = surf_s, surf_e
    389390                   k = surf_usm_h%k(m)
    390 !
    391 !--                Please note, actually, an interpolation of u_0 and v_0
    392 !--                onto the grid center would be required. However, this
    393 !--                would require several data transfers between 2D-grid and
    394 !--                wall type. The effect of this missing interpolation is
    395 !--                negligible. (See also production_e_init).
     391
    396392                   dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
    397393                   dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
     
    405401                DO  m = surf_s, surf_e
    406402                   k = surf_def_h(1)%k(m)
    407 !
    408 !--                Please note, actually, an interpolation of u_0 and v_0
    409 !--                onto the grid center would be required. However, this
    410 !--                would require several data transfers between 2D-grid and
    411 !--                wall type. The effect of this missing interpolation is
    412 !--                negligible. (See also production_e_init).
     403
    413404                   dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
    414405                   dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
    415406
    416407                ENDDO
    417 
    418408             ENDDO
    419409
     
    493483             IF ( .NOT. humidity )  THEN
    494484
    495                 IF ( use_single_reference_value )  THEN
    496 
    497                    IF ( ocean )  THEN
    498 !
    499 !--                   So far in the ocean no special treatment of density flux
    500 !--                   in the bottom and top surface layer
    501                       DO  j = nys, nyn
    502                          DO  k = nzb+1, nzt
    503                             tend(k,j,i) = tend(k,j,i) +                        &
    504                                           kh(k,j,i) * g / rho_reference *      &
    505                                           ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
    506                                           dd2zu(k) *                           &
    507                                     MERGE( 1.0_wp, 0.0_wp,                     &
    508                                              BTEST( wall_flags_0(k,j,i), 0 )   &
    509                                           )
     485                IF ( ocean )  THEN
     486!
     487!--                So far in the ocean no special treatment of density flux
     488!--                in the bottom and top surface layer
     489                   DO  j = nys, nyn
     490                      DO  k = nzb+1, nzt
     491                         tend(k,j,i) = tend(k,j,i) +                           &
     492                                       kh(k,j,i) * g /                         &
     493                              MERGE( rho_reference, prho(k,j,i),               &
     494                                     use_single_reference_value ) *            &
     495                                       ( prho(k+1,j,i) - prho(k-1,j,i) ) *     &
     496                                       dd2zu(k) *                              &
     497                                   MERGE( 1.0_wp, 0.0_wp,                      &
     498                                          BTEST( wall_flags_0(k,j,i), 30 )     &
     499                                        )                            *         &
     500                                   MERGE( 1.0_wp, 0.0_wp,                      &
     501                                          BTEST( wall_flags_0(k,j,i), 9 )      &
     502                                        ) 
     503                      ENDDO
     504!
     505!--                   Treatment of near-surface grid points, at up- and down-
     506!--                   ward facing surfaces
     507                      IF ( use_surface_fluxes )  THEN
     508                         DO  l = 0, 1
     509                            surf_s = surf_def_h(l)%start_index(j,i)
     510                            surf_e = surf_def_h(l)%end_index(j,i)
     511                            DO  m = surf_s, surf_e
     512                               k = surf_def_h(l)%k(m)
     513                               tend(k,j,i) = tend(k,j,i) + g /                 &
     514                                         MERGE( rho_reference, prho(k,j,i),    &
     515                                                use_single_reference_value ) * &
     516                                         drho_air_zw(k-1) *                    &
     517                                         surf_def_h(l)%shf(m)
     518                            ENDDO
    510519                         ENDDO
     520
     521                      ENDIF
     522
     523                      IF ( use_top_fluxes )  THEN
     524                         surf_s = surf_def_h(2)%start_index(j,i)
     525                         surf_e = surf_def_h(2)%end_index(j,i)
     526                         DO  m = surf_s, surf_e
     527                            k = surf_def_h(2)%k(m)
     528                            tend(k,j,i) = tend(k,j,i) + g /                    &
     529                                         MERGE( rho_reference, prho(k,j,i),    &
     530                                                use_single_reference_value ) * &
     531                                         drho_air_zw(k) *                      &
     532                                         surf_def_h(2)%shf(m)
     533                         ENDDO
     534                      ENDIF
     535
     536                   ENDDO
     537
     538                ELSE
     539
     540                   DO  j = nys, nyn
     541                      DO  k = nzb+1, nzt
     542!
     543!--                      Flag 9 is used to mask top fluxes, flag 30 to mask
     544!--                      surface fluxes
     545                         tend(k,j,i) = tend(k,j,i) -                           &
     546                                       kh(k,j,i) * g /                         &
     547                                   MERGE( pt_reference, pt(k,j,i),             &
     548                                           use_single_reference_value ) *      &
     549                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) *         &
     550                                       dd2zu(k)                      *         &
     551                                   MERGE( 1.0_wp, 0.0_wp,                      &
     552                                          BTEST( wall_flags_0(k,j,i), 30 )     &
     553                                        )                            *         &
     554                                   MERGE( 1.0_wp, 0.0_wp,                      &
     555                                          BTEST( wall_flags_0(k,j,i), 9 )      &
     556                                        ) 
    511557                      ENDDO
    512558
    513                    ELSE
    514 
    515                       DO  j = nys, nyn
    516                          DO  k = nzb+1, nzt
    517 !
    518 !--                         Flag 9 is used to mask top fluxes, flag 30 to mask
    519 !--                         surface fluxes
    520                             tend(k,j,i) = tend(k,j,i) -                        &
    521                                           kh(k,j,i) * g / pt_reference  *      &
    522                                           ( pt(k+1,j,i) - pt(k-1,j,i) ) *      &
    523                                           dd2zu(k)                      *      &
    524                                       MERGE( 1.0_wp, 0.0_wp,                   &
    525                                              BTEST( wall_flags_0(k,j,i), 30 )  &
    526                                            )                            *      &
    527                                       MERGE( 1.0_wp, 0.0_wp,                   &
    528                                              BTEST( wall_flags_0(k,j,i), 9 )   &
    529                                            ) 
     559                      IF ( use_surface_fluxes )  THEN
     560!
     561!--                      Default surfaces, up- and downward-facing
     562                         DO  l = 0, 1
     563                            surf_s = surf_def_h(l)%start_index(j,i)
     564                            surf_e = surf_def_h(l)%end_index(j,i)
     565                            DO  m = surf_s, surf_e
     566                               k = surf_def_h(l)%k(m)
     567                               tend(k,j,i) = tend(k,j,i) + g /                 &
     568                                    MERGE( pt_reference, pt(k,j,i),            &
     569                                           use_single_reference_value )        &
     570                                                      * drho_air_zw(k-1)       &
     571                                                      * surf_def_h(l)%shf(m)   
     572                            ENDDO     
    530573                         ENDDO
    531 
    532                          IF ( use_surface_fluxes )  THEN
    533 !
    534 !--                         Default surfaces, up- and downward-facing
    535                             DO  l = 0, 1
    536                                surf_s = surf_def_h(l)%start_index(j,i)
    537                                surf_e = surf_def_h(l)%end_index(j,i)
    538                                DO  m = surf_s, surf_e
    539                                   k = surf_def_h(l)%k(m)
    540                                   tend(k,j,i) = tend(k,j,i) + g / pt_reference &
    541                                                          * drho_air_zw(k-1)    &
    542                                                          * surf_def_h(l)%shf(m)   
    543                                ENDDO   
    544                             ENDDO     
    545 !
    546 !--                         Natural surfaces
    547                             surf_s = surf_lsm_h%start_index(j,i)
    548                             surf_e = surf_lsm_h%end_index(j,i)
    549                             DO  m = surf_s, surf_e
    550                                k = surf_lsm_h%k(m)
    551                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
    552                                                          * drho_air_zw(k-1)    &
    553                                                          * surf_lsm_h%shf(m)   
    554                             ENDDO
    555 !
    556 !--                         Urban surfaces
    557                             surf_s = surf_usm_h%start_index(j,i)
    558                             surf_e = surf_usm_h%end_index(j,i)
    559                             DO  m = surf_s, surf_e
    560                                k = surf_usm_h%k(m)
    561                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
    562                                                          * drho_air_zw(k-1)    &
    563                                                          * surf_usm_h%shf(m)   
    564                             ENDDO                         
    565                          ENDIF
    566 
    567                          IF ( use_top_fluxes )  THEN
    568                             surf_s = surf_def_h(2)%start_index(j,i)
    569                             surf_e = surf_def_h(2)%end_index(j,i)
    570                             DO  m = surf_s, surf_e
    571                                k = surf_def_h(2)%k(m)
    572                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
    573                                                          * drho_air_zw(k-1)    &
    574                                                          * surf_def_h(2)%shf(m)
    575                             ENDDO
    576                          ENDIF
    577                       ENDDO
    578 
    579                    ENDIF
    580 
    581                 ELSE
    582 
    583                    IF ( ocean )  THEN
    584 !
    585 !--                   So far in the ocean no special treatment of density flux
    586 !--                   in the bottom and top surface layer
    587                       DO  j = nys, nyn
    588                          DO  k = nzb+1, nzt
    589                             tend(k,j,i) = tend(k,j,i) +                        &
    590                                           kh(k,j,i) * g / prho(k,j,i) *        &
    591                                           ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
    592                                           dd2zu(k)                           * &
    593                                      MERGE( 1.0_wp, 0.0_wp,                    &
    594                                              BTEST( wall_flags_0(k,j,i), 0 )   &
    595                                           )
     574!
     575!--                      Natural surfaces
     576                         surf_s = surf_lsm_h%start_index(j,i)
     577                         surf_e = surf_lsm_h%end_index(j,i)
     578                         DO  m = surf_s, surf_e
     579                            k = surf_lsm_h%k(m)
     580                            tend(k,j,i) = tend(k,j,i) + g /                    &
     581                                    MERGE( pt_reference, pt(k,j,i),            &
     582                                           use_single_reference_value )        &
     583                                                      * drho_air_zw(k-1)       &
     584                                                      * surf_lsm_h%shf(m)   
    596585                         ENDDO
    597                       ENDDO
    598 
    599                    ELSE
    600 
    601                       DO  j = nys, nyn
    602                          DO  k = nzb+1, nzt
    603 !
    604 !--                         Flag 9 is used to mask top fluxes, flag 30 to mask
    605 !--                         surface fluxes
    606                             tend(k,j,i) = tend(k,j,i) -                        &
    607                                           kh(k,j,i) * g / pt(k,j,i) *          &
    608                                           ( pt(k+1,j,i) - pt(k-1,j,i) ) *      &
    609                                           dd2zu(k)                      *      &
    610                                       MERGE( 1.0_wp, 0.0_wp,                   &
    611                                              BTEST( wall_flags_0(k,j,i), 30 )  &
    612                                            )                            *      &
    613                                       MERGE( 1.0_wp, 0.0_wp,                   &
    614                                              BTEST( wall_flags_0(k,j,i), 9 )   &
    615                                            )
     586!
     587!--                      Urban surfaces
     588                         surf_s = surf_usm_h%start_index(j,i)
     589                         surf_e = surf_usm_h%end_index(j,i)
     590                         DO  m = surf_s, surf_e
     591                            k = surf_usm_h%k(m)
     592                            tend(k,j,i) = tend(k,j,i) + g /                    &
     593                                    MERGE( pt_reference, pt(k,j,i),            &
     594                                           use_single_reference_value )        &
     595                                                      * drho_air_zw(k-1)       &
     596                                                      * surf_usm_h%shf(m)   
     597                         ENDDO                         
     598                      ENDIF
     599
     600                      IF ( use_top_fluxes )  THEN
     601                         surf_s = surf_def_h(2)%start_index(j,i)
     602                         surf_e = surf_def_h(2)%end_index(j,i)
     603                         DO  m = surf_s, surf_e
     604                            k = surf_def_h(2)%k(m)
     605                            tend(k,j,i) = tend(k,j,i) + g /                    &
     606                                    MERGE( pt_reference, pt(k,j,i),            &
     607                                           use_single_reference_value )        &
     608                                                      * drho_air_zw(k)         &
     609                                                      * surf_def_h(2)%shf(m)
    616610                         ENDDO
    617 
    618                          IF ( use_surface_fluxes )  THEN
    619 !
    620 !--                         Default surfaces, up- and downwrd-facing
    621                             DO  l = 0, 1
    622                                surf_s = surf_def_h(l)%start_index(j,i)
    623                                surf_e = surf_def_h(l)%end_index(j,i)
    624                                DO  m = surf_s, surf_e
    625                                   k = surf_def_h(l)%k(m)
    626                                   tend(k,j,i) = tend(k,j,i) + g / pt_reference &
    627                                                          * drho_air_zw(k-1)    &
    628                                                          * surf_def_h(l)%shf(m)   
    629                                ENDDO 
    630                             ENDDO
    631 !
    632 !--                         Natural surfaces
    633                             surf_s = surf_lsm_h%start_index(j,i)
    634                             surf_e = surf_lsm_h%end_index(j,i)
    635                             DO  m = surf_s, surf_e
    636                                k = surf_lsm_h%k(m)
    637                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
    638                                                          * drho_air_zw(k-1)    &
    639                                                          * surf_lsm_h%shf(m)   
    640                             ENDDO 
    641 !
    642 !--                         Urban surfaces
    643                             surf_s = surf_usm_h%start_index(j,i)
    644                             surf_e = surf_usm_h%end_index(j,i)
    645                             DO  m = surf_s, surf_e
    646                                k = surf_usm_h%k(m)
    647                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
    648                                                          * drho_air_zw(k-1)    &
    649                                                          * surf_usm_h%shf(m)   
    650                             ENDDO 
    651                          ENDIF
    652 
    653                          IF ( use_top_fluxes )  THEN
    654                             surf_s = surf_def_h(2)%start_index(j,i)
    655                             surf_e = surf_def_h(2)%end_index(j,i)
    656                             DO  m = surf_s, surf_e
    657                                k = surf_def_h(2)%k(m)
    658                                tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
    659                                                          * drho_air_zw(k-1)    &
    660                                                          * surf_def_h(2)%shf(m)
    661                             ENDDO
    662                          ENDIF
    663                       ENDDO
    664 
    665                    ENDIF
     611                      ENDIF
     612                   ENDDO
    666613
    667614                ENDIF
     
    679626                         k2 = 0.61_wp * pt(k,j,i)
    680627                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
    681                                          g / vpt(k,j,i) *                      &
     628                                         g /                                   &
     629                                    MERGE( vpt_reference, vpt(k,j,i),          &
     630                                           use_single_reference_value ) *      &
    682631                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
    683632                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
     
    704653                         ENDIF
    705654                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
    706                                          g / vpt(k,j,i) *                      &
     655                                         g /                                   &
     656                                    MERGE( vpt_reference, vpt(k,j,i),          &
     657                                           use_single_reference_value ) *      &
    707658                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
    708659                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
     
    718669                         k2 = 0.61_wp * pt(k,j,i)
    719670                         tend(k,j,i) = tend(k,j,i) -                           &
    720                                        kh(k,j,i) * g / vpt(k,j,i) *            &
     671                                       kh(k,j,i) * g /                         &
     672                                    MERGE( vpt_reference, vpt(k,j,i),          &
     673                                           use_single_reference_value ) *      &
    721674                                       ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +   &
    722675                                         k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -   &
     
    739692                   DO  j = nys, nyn
    740693!
    741 !--                   Treat horizontal default surfaces, up- and downward-facing
     694!--                   Treat horizontal default surfaces
    742695                      DO  l = 0, 1
    743696                         surf_s = surf_def_h(l)%start_index(j,i)
     
    768721                            ENDIF
    769722
    770                             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *       &
     723                            tend(k,j,i) = tend(k,j,i) + g /                    &
     724                                    MERGE( vpt_reference, vpt(k,j,i),          &
     725                                           use_single_reference_value ) *      &
    771726                                               ( k1 * surf_def_h(l)%shf(m) +   &
    772                                                  k2 * surf_def_h(l)%qsws(m)       &
     727                                                 k2 * surf_def_h(l)%qsws(m)    &
    773728                                               ) * drho_air_zw(k-1)
    774729                         ENDDO
     
    803758                         ENDIF
    804759
    805                          tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
     760                         tend(k,j,i) = tend(k,j,i) + g /                       &
     761                                    MERGE( vpt_reference, vpt(k,j,i),          &
     762                                           use_single_reference_value ) *      &
    806763                                               ( k1 * surf_lsm_h%shf(m) +      &
    807764                                                 k2 * surf_lsm_h%qsws(m)       &
     
    837794                         ENDIF
    838795
    839                          tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
     796                         tend(k,j,i) = tend(k,j,i) + g /                       &
     797                                    MERGE( vpt_reference, vpt(k,j,i),          &
     798                                           use_single_reference_value ) *      &
    840799                                               ( k1 * surf_usm_h%shf(m) +      &
    841800                                                 k2 * surf_usm_h%qsws(m)       &
     
    878837                         ENDIF
    879838
    880                          tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
     839                         tend(k,j,i) = tend(k,j,i) + g /                       &
     840                                    MERGE( vpt_reference, vpt(k,j,i),          &
     841                                           use_single_reference_value ) *      &
    881842                                               ( k1 * surf_def_h(2)%shf(m) +   &
    882                                                  k2 * surf_def_h(2)%qsws(m)       &
    883                                                ) * drho_air_zw(k-1)
     843                                                 k2 * surf_def_h(2)%qsws(m)    &
     844                                               ) * drho_air_zw(k)
    884845
    885846                      ENDDO
     
    916877                  humidity, kappa, neutral, ocean, pt_reference,               &
    917878                  rho_reference, use_single_reference_value,                   &
    918                   use_surface_fluxes, use_top_fluxes
     879                  use_surface_fluxes, use_top_fluxes, vpt_reference
    919880
    920881       USE grid_variables,                                                     &
     
    10651026!--             -1.0 for right-facing wall, 1.0 for left-facing wall
    10661027                sign_dir = MERGE( 1.0_wp, -1.0_wp,                             &
    1067                                   BTEST( wall_flags_0(k,j,i-1), 0 ) ) 
     1028                                  BTEST( wall_flags_0(k,j,i-1), 0 ) )
    10681029                dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    10691030                dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
     
    11291090          DO  m = surf_s, surf_e
    11301091             k = surf_lsm_h%k(m)
    1131 !
    1132 !--          Please note, actually, an interpolation of u_0 and v_0
    1133 !--          onto the grid center would be required. However, this
    1134 !--          would require several data transfers between 2D-grid and
    1135 !--          wall type. The effect of this missing interpolation is
    1136 !--          negligible. (See also production_e_init).
     1092
    11371093             dudz(k)     = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
    11381094             dvdz(k)     = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
     
    11441100          DO  m = surf_s, surf_e
    11451101             k = surf_usm_h%k(m)
    1146 !
    1147 !--          Please note, actually, an interpolation of u_0 and v_0
    1148 !--          onto the grid center would be required. However, this
    1149 !--          would require several data transfers between 2D-grid and
    1150 !--          wall type. The effect of this missing interpolation is
    1151 !--          negligible. (See also production_e_init).
     1102
    11521103             dudz(k)     = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
    11531104             dvdz(k)     = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
     
    11601111          DO  m = surf_s, surf_e
    11611112             k = surf_def_h(1)%k(m)
    1162 !
    1163 !--          Please note, actually, an interpolation of u_0 and v_0
    1164 !--          onto the grid center would be required. However, this
    1165 !--          would require several data transfers between 2D-grid and
    1166 !--          wall type. The effect of this missing interpolation is
    1167 !--          negligible. (See also production_e_init).
     1113
    11681114             dudz(k)     = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
    11691115             dvdz(k)     = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
     
    12321178          IF ( .NOT. humidity )  THEN
    12331179
    1234              IF ( use_single_reference_value )  THEN
    1235 
    1236                 IF ( ocean )  THEN
    1237 !
    1238 !--                So far in the ocean no special treatment of density flux in
    1239 !--                the bottom and top surface layer
    1240                    DO  k = nzb+1, nzt
     1180             IF ( ocean )  THEN
     1181!
     1182!--             So far in the ocean no special treatment of density flux in
     1183!--             the bottom and top surface layer
     1184                DO  k = nzb+1, nzt
    12411185 
    1242                       tend(k,j,i) = tend(k,j,i) +                              &
    1243                                     kh(k,j,i) * g / rho_reference *            &
    1244                                     ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
    1245                                     dd2zu(k) *                                 &
    1246                                       MERGE( 1.0_wp, 0.0_wp,                   &
    1247                                              BTEST( wall_flags_0(k,j,i), 0 )   &
    1248                                            )                   
     1186                   tend(k,j,i) = tend(k,j,i) +                                 &
     1187                                 kh(k,j,i) * g /                               &
     1188                                 MERGE( rho_reference, prho(k,j,i),            &
     1189                                        use_single_reference_value ) *         &
     1190                                 ( prho(k+1,j,i) - prho(k-1,j,i) ) *           &
     1191                                 dd2zu(k) *                                    &
     1192                                   MERGE( 1.0_wp, 0.0_wp,                      &
     1193                                          BTEST( wall_flags_0(k,j,i), 30 )     &
     1194                                        ) *                                    &
     1195                                   MERGE( 1.0_wp, 0.0_wp,                      &
     1196                                          BTEST( wall_flags_0(k,j,i), 9 )      &
     1197                                        )           
     1198                ENDDO
     1199
     1200                IF ( use_surface_fluxes )  THEN
     1201!
     1202!--                Default surfaces, up- and downward-facing
     1203                   DO  l = 0, 1
     1204                      surf_s = surf_def_h(l)%start_index(j,i)
     1205                      surf_e = surf_def_h(l)%end_index(j,i)
     1206                      DO  m = surf_s, surf_e
     1207                         k = surf_def_h(l)%k(m)
     1208                         tend(k,j,i) = tend(k,j,i) + g /                       &
     1209                                   MERGE( rho_reference, prho(k,j,i),          &
     1210                                          use_single_reference_value ) *       &
     1211                                   drho_air_zw(k-1) *                          &
     1212                                   surf_def_h(l)%shf(m)
     1213                      ENDDO
    12491214                   ENDDO
    12501215
    1251                 ELSE
    1252 
    1253                    DO  k = nzb+1, nzt
    1254 !
    1255 !--                   Flag 9 is used to mask top fluxes, flag 30 to mask
    1256 !--                   surface fluxes
    1257                       tend(k,j,i) = tend(k,j,i) -                              &
    1258                                     kh(k,j,i) * g / pt_reference *             &
    1259                                     ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
    1260                                       MERGE( 1.0_wp, 0.0_wp,                   &
    1261                                              BTEST( wall_flags_0(k,j,i), 30 )  &
    1262                                            ) *                                 &
    1263                                       MERGE( 1.0_wp, 0.0_wp,                   &
    1264                                              BTEST( wall_flags_0(k,j,i), 9 )   &
    1265                                            )
    1266 
     1216                ENDIF
     1217
     1218                IF ( use_top_fluxes )  THEN
     1219                   surf_s = surf_def_h(2)%start_index(j,i)
     1220                   surf_e = surf_def_h(2)%end_index(j,i)
     1221                   DO  m = surf_s, surf_e
     1222                      k = surf_def_h(2)%k(m)
     1223                      tend(k,j,i) = tend(k,j,i) + g /                          &
     1224                                   MERGE( rho_reference, prho(k,j,i),          &
     1225                                          use_single_reference_value ) *       &
     1226                                   drho_air_zw(k) *                            &
     1227                                   surf_def_h(2)%shf(m)
    12671228                   ENDDO
    1268 
    1269                    IF ( use_surface_fluxes )  THEN
    1270 !
    1271 !--                   Default surfaces, up- and downward-facing
    1272                       DO  l = 0, 1
    1273                          surf_s = surf_def_h(l)%start_index(j,i)
    1274                          surf_e = surf_def_h(l)%end_index(j,i)
    1275                          DO  m = surf_s, surf_e
    1276                             k = surf_def_h(l)%k(m)
    1277                             tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
    1278                                                            drho_air_zw(k-1) *  &
    1279                                                            surf_def_h(l)%shf(m)
    1280                          ENDDO
     1229                ENDIF
     1230
     1231             ELSE
     1232
     1233                DO  k = nzb+1, nzt
     1234!
     1235!--                Flag 9 is used to mask top fluxes, flag 30 to mask
     1236!--                surface fluxes
     1237                   tend(k,j,i) = tend(k,j,i) -                                 &
     1238                                 kh(k,j,i) * g /                               &
     1239                                   MERGE( pt_reference, pt(k,j,i),             &
     1240                                          use_single_reference_value ) *       &
     1241                                 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) *    &
     1242                                   MERGE( 1.0_wp, 0.0_wp,                      &
     1243                                          BTEST( wall_flags_0(k,j,i), 30 )     &
     1244                                        ) *                                    &
     1245                                   MERGE( 1.0_wp, 0.0_wp,                      &
     1246                                          BTEST( wall_flags_0(k,j,i), 9 )      &
     1247                                        )
     1248
     1249                ENDDO
     1250
     1251                IF ( use_surface_fluxes )  THEN
     1252!
     1253!--                Default surfaces, up- and downward-facing
     1254                   DO  l = 0, 1
     1255                      surf_s = surf_def_h(l)%start_index(j,i)
     1256                      surf_e = surf_def_h(l)%end_index(j,i)
     1257                      DO  m = surf_s, surf_e
     1258                         k = surf_def_h(l)%k(m)
     1259                         tend(k,j,i) = tend(k,j,i) + g /                       &
     1260                                   MERGE( pt_reference, pt(k,j,i),             &
     1261                                          use_single_reference_value ) *       &
     1262                                   drho_air_zw(k-1) *                          &
     1263                                   surf_def_h(l)%shf(m)
    12811264                      ENDDO
    1282 !
    1283 !--                   Natural surfaces
    1284                       surf_s = surf_lsm_h%start_index(j,i)
    1285                       surf_e = surf_lsm_h%end_index(j,i)
    1286                       DO  m = surf_s, surf_e
    1287                          k = surf_lsm_h%k(m)
    1288                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
    1289                                                            drho_air_zw(k-1) *  &
    1290                                                            surf_lsm_h%shf(m)
    1291                       ENDDO
    1292 !
    1293 !--                   Urban surfaces
    1294                       surf_s = surf_usm_h%start_index(j,i)
    1295                       surf_e = surf_usm_h%end_index(j,i)
    1296                       DO  m = surf_s, surf_e
    1297                          k = surf_usm_h%k(m)
    1298                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
    1299                                                            drho_air_zw(k-1) *  &
    1300                                                            surf_usm_h%shf(m)
    1301                       ENDDO
    1302                    ENDIF
    1303 
    1304                    IF ( use_top_fluxes )  THEN
    1305                       surf_s = surf_def_h(2)%start_index(j,i)
    1306                       surf_e = surf_def_h(2)%end_index(j,i)
    1307                       DO  m = surf_s, surf_e
    1308                          k = surf_def_h(2)%k(m)
    1309                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
    1310                                                            drho_air_zw(k-1) *  &
    1311                                                            surf_def_h(2)%shf(m)
    1312                       ENDDO
    1313                    ENDIF
    1314 
     1265                   ENDDO
     1266!
     1267!--                Natural surfaces
     1268                   surf_s = surf_lsm_h%start_index(j,i)
     1269                   surf_e = surf_lsm_h%end_index(j,i)
     1270                   DO  m = surf_s, surf_e
     1271                      k = surf_lsm_h%k(m)
     1272                      tend(k,j,i) = tend(k,j,i) + g /                          &
     1273                                   MERGE( pt_reference, pt(k,j,i),             &
     1274                                          use_single_reference_value ) *       &
     1275                                   drho_air_zw(k-1) *                          &
     1276                                   surf_lsm_h%shf(m)
     1277                   ENDDO
     1278!
     1279!--                Urban surfaces
     1280                   surf_s = surf_usm_h%start_index(j,i)
     1281                   surf_e = surf_usm_h%end_index(j,i)
     1282                   DO  m = surf_s, surf_e
     1283                      k = surf_usm_h%k(m)
     1284                      tend(k,j,i) = tend(k,j,i) + g /                          &
     1285                                   MERGE( pt_reference, pt(k,j,i),             &
     1286                                          use_single_reference_value ) *       &
     1287                                   drho_air_zw(k-1) *                          &
     1288                                   surf_usm_h%shf(m)
     1289                   ENDDO
    13151290                ENDIF
    13161291
    1317              ELSE
    1318 
    1319                 IF ( ocean )  THEN
    1320 !
    1321 !--                So far in the ocean no special treatment of density flux in
    1322 !--                the bottom and top surface layer
    1323                    DO  k = nzb+1, nzt
    1324                       tend(k,j,i) = tend(k,j,i) +                              &
    1325                                     kh(k,j,i) * g / prho(k,j,i)       *        &
    1326                                     ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
    1327                                     dd2zu(k) *                                 &
    1328                                       MERGE( 1.0_wp, 0.0_wp,                   &
    1329                                              BTEST( wall_flags_0(k,j,i), 0 )   &
    1330                                            )           
     1292                IF ( use_top_fluxes )  THEN
     1293                   surf_s = surf_def_h(2)%start_index(j,i)
     1294                   surf_e = surf_def_h(2)%end_index(j,i)
     1295                   DO  m = surf_s, surf_e
     1296                      k = surf_def_h(2)%k(m)
     1297                      tend(k,j,i) = tend(k,j,i) + g /                          &
     1298                                   MERGE( pt_reference, pt(k,j,i),             &
     1299                                          use_single_reference_value ) *       &
     1300                                   drho_air_zw(k) *                            &
     1301                                   surf_def_h(2)%shf(m)
    13311302                   ENDDO
    1332 
    1333                 ELSE
    1334 
    1335                    DO  k = nzb+1, nzt
    1336 !
    1337 !--                   Flag 9 is used to mask top fluxes, flag 30 to mask
    1338 !--                   surface fluxes
    1339                       tend(k,j,i) = tend(k,j,i) -                              &
    1340                                     kh(k,j,i) * g / pt(k,j,i) *                &
    1341                                     ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
    1342                                       MERGE( 1.0_wp, 0.0_wp,                   &
    1343                                              BTEST( wall_flags_0(k,j,i), 30 )  &
    1344                                            ) *                                 &
    1345                                       MERGE( 1.0_wp, 0.0_wp,                   &
    1346                                              BTEST( wall_flags_0(k,j,i), 9 )   &
    1347                                            )
    1348                    ENDDO
    1349 
    1350                    IF ( use_surface_fluxes )  THEN
    1351 !
    1352 !--                   Default surfaces, up- and downward-facing
    1353                       DO  l = 0, 1
    1354                          surf_s = surf_def_h(l)%start_index(j,i)
    1355                          surf_e = surf_def_h(l)%end_index(j,i)
    1356                          DO  m = surf_s, surf_e
    1357                             k = surf_def_h(l)%k(m)
    1358                             tend(k,j,i) = tend(k,j,i) + g / pt_reference       &
    1359                                                    * drho_air_zw(k-1)          &
    1360                                                    * surf_def_h(l)%shf(m)
    1361                          ENDDO 
    1362                       ENDDO
    1363 !
    1364 !--                   Natural surfaces
    1365                       surf_s = surf_lsm_h%start_index(j,i)
    1366                       surf_e = surf_lsm_h%end_index(j,i)
    1367                       DO  m = surf_s, surf_e
    1368                          k = surf_lsm_h%k(m)
    1369                          tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
    1370                                                      * drho_air_zw(k-1)        &
    1371                                                      * surf_lsm_h%shf(m)
    1372                       ENDDO 
    1373 !
    1374 !--                   Urban surfaces
    1375                       surf_s = surf_usm_h%start_index(j,i)
    1376                       surf_e = surf_usm_h%end_index(j,i)
    1377                       DO  m = surf_s, surf_e
    1378                          k = surf_usm_h%k(m)
    1379                          tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
    1380                                                      * drho_air_zw(k-1)        &
    1381                                                      * surf_usm_h%shf(m)
    1382                       ENDDO
    1383                    ENDIF
    1384 
    1385                    IF ( use_top_fluxes )  THEN
    1386                       surf_s = surf_def_h(2)%start_index(j,i)
    1387                       surf_e = surf_def_h(2)%end_index(j,i)
    1388                       DO  m = surf_s, surf_e
    1389                          k = surf_def_h(2)%k(m)
    1390                          tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
    1391                                                            drho_air_zw(k-1) *  &
    1392                                                            surf_def_h(2)%shf(m)
    1393                       ENDDO
    1394                    ENDIF
    1395 
    13961303                ENDIF
    13971304
     
    14061313                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    14071314                   k2 = 0.61_wp * pt(k,j,i)
    1408                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
     1315                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
     1316                                   MERGE( vpt_reference, vpt(k,j,i),           &
     1317                                          use_single_reference_value ) *       &
    14091318                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
    14101319                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
     
    14311340                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    14321341                   ENDIF
    1433                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
     1342                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
     1343                                   MERGE( vpt_reference, vpt(k,j,i),           &
     1344                                          use_single_reference_value ) *       &
    14341345                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
    14351346                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
     
    14441355                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
    14451356                   k2 = 0.61_wp * pt(k,j,i)
    1446                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *    &
     1357                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                 &
     1358                                   MERGE( vpt_reference, vpt(k,j,i),           &
     1359                                          use_single_reference_value ) *       &
    14471360                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +      &
    14481361                                       k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -      &
     
    14891402                      ENDIF
    14901403
    1491                       tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *             &
     1404                      tend(k,j,i) = tend(k,j,i) + g /                          &
     1405                                   MERGE( vpt_reference, vpt(k,j,i),           &
     1406                                          use_single_reference_value ) *       &
    14921407                                         ( k1 * surf_def_h(l)%shf(m) +         &
    14931408                                           k2 * surf_def_h(l)%qsws(m)          &
     
    15241439                   ENDIF
    15251440
    1526                    tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
     1441                   tend(k,j,i) = tend(k,j,i) + g /                             &
     1442                                   MERGE( vpt_reference, vpt(k,j,i),           &
     1443                                          use_single_reference_value ) *       &
    15271444                                            ( k1 * surf_lsm_h%shf(m) +         &
    15281445                                              k2 * surf_lsm_h%qsws(m)          &
     
    15581475                   ENDIF
    15591476
    1560                    tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
     1477                   tend(k,j,i) = tend(k,j,i) + g /                             &
     1478                                   MERGE( vpt_reference, vpt(k,j,i),           &
     1479                                          use_single_reference_value ) *       &
    15611480                                            ( k1 * surf_usm_h%shf(m) +         &
    15621481                                              k2 * surf_usm_h%qsws(m)          &
     
    15961515                   ENDIF
    15971516
    1598                    tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
     1517                   tend(k,j,i) = tend(k,j,i) + g /                             &
     1518                                   MERGE( vpt_reference, vpt(k,j,i),           &
     1519                                          use_single_reference_value ) *       &
    15991520                               ( k1* surf_def_h(2)%shf(m) +                    &
    16001521                                 k2 * surf_def_h(2)%qsws(m)                    &
    1601                                ) * drho_air_zw(k-1)
     1522                               ) * drho_air_zw(k)
    16021523                ENDDO
    16031524
     
    16821603          ENDDO
    16831604!
    1684 !--       Default surfaces, downward-facing
     1605!--       Default surfaces, downward-facing surfaces
    16851606          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
    16861607          DO  m = 1, surf_def_h(1)%ns
     
    16891610             j = surf_def_h(1)%j(m)
    16901611             k = surf_def_h(1)%k(m)
    1691 !
    1692 !--          Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
    1693 !--          and km are not on the same grid. Actually, a further
    1694 !--          interpolation of km onto the u/v-grid is necessary. However, the
    1695 !--          effect of this error is negligible.
    1696 !--          In case of downward-facing surfaces, gradient is calculated
    1697 !--          between u_0 and u(k-1).
    1698              surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *       &
    1699                                         drho_air_zw(k-1) *                     &
    1700                                         ( zu(k+1)    - zu(k-1)    )  /         &
     1612
     1613             surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *    &
     1614                                        drho_air_zw(k-1) *                  &
     1615                                        ( zu(k+1)    - zu(k-1)    )  /      &
    17011616                                        ( km(k,j,i)  + 1.0E-20_wp ) 
    1702              surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *       &
    1703                                         drho_air_zw(k-1) *                     &
    1704                                         ( zu(k+1)    - zu(k-1)    )  /         &
     1617             surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *    &
     1618                                        drho_air_zw(k-1) *                  &
     1619                                        ( zu(k+1)    - zu(k-1)    )  /      &
    17051620                                        ( km(k,j,i)  + 1.0E-20_wp ) 
    17061621
    1707              IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >                  &
    1708                   ABS( u(k+1,j,i)           - u(k-1,j,i) )                     &
     1622             IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >               &
     1623                  ABS( u(k+1,j,i)           - u(k-1,j,i) )                  &
    17091624                )  surf_def_h(1)%u_0(m) = u(k+1,j,i)
    17101625
    1711              IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >                  &
    1712                   ABS( v(k+1,j,i)           - v(k-1,j,i) )                     &
     1626             IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >               &
     1627                  ABS( v(k+1,j,i)           - v(k-1,j,i) )                  &
    17131628                )  surf_def_h(1)%v_0(m) = v(k+1,j,i)
    17141629
  • palm/trunk/SOURCE/surface_mod.f90

    r2478 r2508  
    2525! -----------------
    2626! $Id$
     27! Minor formatting adjustment
     28!
     29! 2478 2017-09-18 13:37:24Z suehring
    2730! Bugfixes in initializing model top
    2831!
     
    795798       ALLOCATE ( surfaces%j(1:surfaces%ns)  )
    796799       ALLOCATE ( surfaces%k(1:surfaces%ns)  )
     800
     801       IF ( .NOT. constant_diffusion )  THEN   
     802          ALLOCATE ( surfaces%u_0(1:surfaces%ns) ) 
     803          ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
     804       ENDIF
    797805!
    798806!--    Vertical momentum fluxes of u and v
     
    15041512!--          Prescribe latent heat flux at the top     
    15051513             IF ( humidity )  THEN
    1506              surf%qsws(num_h) = 0.0_wp
     1514                surf%qsws(num_h) = 0.0_wp
    15071515                IF ( cloud_physics  .AND.  microphysics_morrison ) THEN
    15081516                   surf%ncsws(num_h) = 0.0_wp
Note: See TracChangeset for help on using the changeset viewer.