Ignore:
Timestamp:
Nov 4, 2015 2:47:01 PM (6 years ago)
Author:
maronga
Message:

several bugfixes related to the new surface layer routine and land-surface-radiation interaction

File:
1 edited

Legend:

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

    r1706 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Bugfix: division by zero could occur when calculating rib at low wind speeds
     22! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
     23! neutral = .T.
    2224!
    2325! Former revisions:
     
    9698!>    that this method cannot be used in case of roughness heterogeneity
    9799!>
    98 !> @todo limiting of uv_total might be required in some cases
    99100!> @todo (re)move large_scale_forcing actions
    100101!> @todo check/optimize OpenMP and OpenACC directives
     
    208209          CALL calc_scaling_parameters
    209210
     211          CALL calc_uv_total
     212
    210213          IF ( .NOT. neutral )  THEN
    211214             CALL calc_ol
     
    220223!--    number to calculate the Obukhov length
    221224       ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
     225
     226          CALL calc_uv_total
    222227
    223228          IF ( .NOT. neutral )  THEN
     
    273278!--    Allocate field for storing the horizontal velocity
    274279       ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
     280
     281
     282!
     283!--    In case of runs with neutral statification, set Obukhov length to a
     284!--    large value
     285       IF ( neutral ) ol = 1.0E10_wp
    275286
    276287       IF ( most_method == 'lookup' )  THEN
     
    384395! Description:
    385396! ------------
     397!> Compute the absolute value of the horizontal velocity (relative to the
     398!> surface). This is required by all methods
     399!------------------------------------------------------------------------------!
     400    SUBROUTINE calc_uv_total
     401
     402       IMPLICIT NONE
     403
     404
     405       !$OMP PARALLEL DO PRIVATE( k, z_mo )
     406       !$acc kernels loop
     407       DO  i = nxl, nxr
     408          DO  j = nys, nyn
     409
     410             k   = nzb_s_inner(j,i)
     411             uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
     412                                         - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
     413                              ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
     414                                         - v(k,j,i)   - v(k,j+1,i) ) )**2 )
     415
     416!
     417!--          For too small values of the local wind, MOST does not work. A
     418!--          threshold value is thus set if required
     419!            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
     420
     421          ENDDO
     422       ENDDO
     423
     424!
     425!--    Values of uv_total need to be exchanged at the ghost boundaries
     426       !$acc update host( uv_total )
     427       CALL exchange_horiz_2d( uv_total )
     428       !$acc update device( uv_total )
     429
     430    END SUBROUTINE calc_uv_total
     431
     432
     433!------------------------------------------------------------------------------!
     434! Description:
     435! ------------
    386436!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
    387437!------------------------------------------------------------------------------!
     
    394444
    395445       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
    396                        
     446
    397447       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
    398448                       f_d_ol, & !< Derivative of f
    399449                       ol_l,   & !< Lower bound of L for Newton iteration
    400450                       ol_m,   & !< Previous value of L for Newton iteration
    401                        ol_old, & !< Previous time step value of L 
     451                       ol_old, & !< Previous time step value of L
    402452                       ol_u      !< Upper bound of L for Newton iteration
    403453
    404 !
    405 !--    Compute the absolute value of the horizontal velocity (relative to the
    406 !--    surface). This is required by all methods
    407        !$OMP PARALLEL DO PRIVATE( k, z_mo )
    408        !$acc kernels loop
    409        DO  i = nxl, nxr
    410           DO  j = nys, nyn
    411 
    412              k   = nzb_s_inner(j,i)
    413 
    414              uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
    415                                          - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
    416                               ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
    417                                          - v(k,j,i)   - v(k,j+1,i) ) )**2 )
    418 
    419 !
    420 !--          For too small values of the local wind, MOST does not work. A
    421 !--          threshold value is thus set if required
    422 !            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
    423 
    424           ENDDO
    425        ENDDO
    426 
    427 !
    428 !--    Values of uv_total need to be exchanged at the ghost boundaries
    429        !$acc update host( uv_total )
    430        CALL exchange_horiz_2d( uv_total )
    431        !$acc update device( uv_total )
    432454
    433455       IF ( TRIM( most_method ) /= 'circular' )  THEN
     
    447469                   IF ( humidity )  THEN
    448470                      rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) )      &
    449                            / ( uv_total(j,i)**2 * vpt(k+1,j,i) )
     471                           / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E-20_wp )
    450472                   ELSE
    451473                      rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) )        &
    452                            / ( uv_total(j,i)**2 * pt(k+1,j,i) )
     474                           / ( uv_total(j,i)**2 * pt(k+1,j,i)  + 1.0E-20_wp )
    453475                   ENDIF     
    454476                ELSE
     
    462484                                 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
    463485                                 * pt(k+1,j,i) * qsws(j,i) )   &
    464                                  / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2 )
     486                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
     487                                 + 1.0E-20_wp)
    465488                   ELSE
    466489                      rib(j,i) = - g * z_mo * shf(j,i)                         &
    467                            / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 )
     490                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
     491                           + 1.0E-20_wp )
    468492                   ENDIF
    469493                ENDIF 
Note: See TracChangeset for help on using the changeset viewer.