Ignore:
Timestamp:
Jan 14, 2019 12:49:24 PM (5 years ago)
Author:
maronga
Message:

removed most_methods circular and lookup. added improved version of palm_csd

File:
1 edited

Legend:

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

    r3655 r3668  
    2626! -----------------
    2727! $Id$
     28! Removed methods "circular" and "lookup"
     29!
     30! 3655 2019-01-07 16:51:22Z knoop
    2831! OpenACC port for SPEC
    2932!
     
    228231! ------------
    229232!> Diagnostic computation of vertical fluxes in the constant flux layer from the
    230 !> values of the variables at grid point k=1. Three different methods are
    231 !> available:
    232 !> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
    233 !> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
    234 !>    slower
    235 !> 3) a method using a lookup table which is fast and accurate. Note, however,
    236 !>    that this method cannot be used in case of roughness heterogeneity
     233!> values of the variables at grid point k=1 based on Newton iteration
    237234!>
    238235!> @todo (re)move large_scale_forcing actions
     
    261258               intermediate_timestep_count, intermediate_timestep_count_max,   &
    262259               land_surface, large_scale_forcing, lsf_surf, message_string,    &
    263                most_method, neutral, passive_scalar, pt_surface, q_surface,    &
     260               neutral, passive_scalar, pt_surface, q_surface,                 &
    264261               run_coupled, surface_pressure, simulated_time, terminate_run,   &
    265262               time_since_reference_point, urban_surface,                      &
     
    400397!--    First, calculate the new Obukhov length, then new friction velocity,
    401398!--    followed by the new scaling parameters (th*, q*, etc.), and the new
    402 !--    surface fluxes if required. The old routine ("circular") requires a
    403 !--    different order of calls as the scaling parameters from the previous time
    404 !--    steps are used to calculate the Obukhov length
    405 
    406 !
    407 !--    Depending on setting of most_method use the "old" routine
    408 !--    Note, each routine is called for different surface types.
     399!--    surface fluxes if required. Note, each routine is called for different surface types.
    409400!--    First call for default-type horizontal surfaces, for natural- and
    410401!--    urban-type surfaces. Note, at this place only upward-facing horizontal
    411 !--    surfaces are treted.
    412        IF ( most_method == 'circular' )  THEN
    413 !
    414 !--       Default-type upward-facing horizontal surfaces
    415           IF ( surf_def_h(0)%ns >= 1  )  THEN
    416              surf => surf_def_h(0)
    417              CALL calc_scaling_parameters
    418              CALL calc_uvw_abs
    419              IF ( .NOT. neutral )  CALL calc_ol
    420              CALL calc_us
    421              CALL calc_surface_fluxes
    422              
    423              IF ( do_output_at_2m )  THEN
    424                 CALL calc_pt_near_surface ( '2m' )
    425              ENDIF
    426           ENDIF
    427 !
    428 !--       Natural-type horizontal surfaces
    429           IF ( surf_lsm_h%ns >= 1  )  THEN
    430              surf => surf_lsm_h
    431              CALL calc_scaling_parameters
    432              CALL calc_uvw_abs
    433              IF ( .NOT. neutral )  CALL calc_ol
    434              CALL calc_us
    435              CALL calc_surface_fluxes
    436              IF ( do_output_at_2m )  THEN
    437                 CALL calc_pt_near_surface ( '2m' )
    438              ENDIF
    439           ENDIF
    440 !
    441 !--       Urban-type horizontal surfaces
    442           IF ( surf_usm_h%ns >= 1  )  THEN
    443              surf => surf_usm_h
    444              CALL calc_scaling_parameters
    445              CALL calc_uvw_abs
    446              IF ( .NOT. neutral )  CALL calc_ol
    447              CALL calc_us
    448              CALL calc_surface_fluxes
    449              IF ( do_output_at_2m )  THEN
    450                 CALL calc_pt_near_surface ( '2m' )
    451              ENDIF
    452              IF ( indoor_model )  THEN
    453                 CALL calc_pt_near_surface ( '10cm' )
    454              ENDIF
    455           ENDIF
    456 !
    457 !--    Use either Newton iteration or a lookup table for the bulk Richardson
    458 !--    number to calculate the Obukhov length
    459        ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
    460 !
    461 !--       Default-type upward-facing horizontal surfaces
    462           IF ( surf_def_h(0)%ns >= 1  )  THEN
    463              surf => surf_def_h(0)
    464              CALL calc_uvw_abs
    465              IF ( .NOT. neutral )  CALL calc_ol
    466              CALL calc_us
    467              CALL calc_scaling_parameters
    468              CALL calc_surface_fluxes
    469              IF ( do_output_at_2m )  THEN
    470                 CALL calc_pt_near_surface ( '2m' )
    471              ENDIF
    472           ENDIF
    473 !
    474 !--       Natural-type horizontal surfaces
    475           IF ( surf_lsm_h%ns >= 1  )  THEN
    476              surf => surf_lsm_h
    477              CALL calc_uvw_abs
    478              IF ( .NOT. neutral )  CALL calc_ol
    479              CALL calc_us
    480              CALL calc_scaling_parameters
    481              CALL calc_surface_fluxes
    482              IF ( do_output_at_2m )  THEN
    483                 CALL calc_pt_near_surface ( '2m' )
    484              ENDIF
    485           ENDIF
    486 !
    487 !--       Urban-type horizontal surfaces
    488           IF ( surf_usm_h%ns >= 1  )  THEN
    489              surf => surf_usm_h
    490              CALL calc_uvw_abs
    491              IF ( .NOT. neutral )  CALL calc_ol
    492              CALL calc_us
    493              CALL calc_scaling_parameters
    494              CALL calc_surface_fluxes
    495              IF ( do_output_at_2m )  THEN
    496                 CALL calc_pt_near_surface ( '2m' )
    497              ENDIF
    498              IF ( indoor_model )  THEN
    499                 CALL calc_pt_near_surface ( '10cm' )
    500              ENDIF
    501           ENDIF
    502 
     402!--    surfaces are treated.
     403
     404!
     405!--    Default-type upward-facing horizontal surfaces
     406       IF ( surf_def_h(0)%ns >= 1  )  THEN
     407          surf => surf_def_h(0)
     408          CALL calc_uvw_abs
     409          IF ( .NOT. neutral )  CALL calc_ol
     410          CALL calc_us
     411          CALL calc_scaling_parameters
     412          CALL calc_surface_fluxes
     413          IF ( do_output_at_2m )  THEN
     414             CALL calc_pt_near_surface ( '2m' )
     415          ENDIF
    503416       ENDIF
     417!
     418!--    Natural-type horizontal surfaces
     419       IF ( surf_lsm_h%ns >= 1  )  THEN
     420          surf => surf_lsm_h
     421          CALL calc_uvw_abs
     422          IF ( .NOT. neutral )  CALL calc_ol
     423          CALL calc_us
     424          CALL calc_scaling_parameters
     425          CALL calc_surface_fluxes
     426          IF ( do_output_at_2m )  THEN
     427             CALL calc_pt_near_surface ( '2m' )
     428          ENDIF
     429       ENDIF
     430!
     431!--    Urban-type horizontal surfaces
     432       IF ( surf_usm_h%ns >= 1  )  THEN
     433          surf => surf_usm_h
     434          CALL calc_uvw_abs
     435          IF ( .NOT. neutral )  CALL calc_ol
     436          CALL calc_us
     437          CALL calc_scaling_parameters
     438          CALL calc_surface_fluxes
     439          IF ( do_output_at_2m )  THEN
     440             CALL calc_pt_near_surface ( '2m' )
     441          ENDIF
     442          IF ( indoor_model )  THEN
     443             CALL calc_pt_near_surface ( '10cm' )
     444          ENDIF
     445       ENDIF
     446
    504447!
    505448!--    Treat downward-facing horizontal surfaces. Note, so far, these are
     
    549492!--    flux in land-surface model in case of aero_resist_kray is not true.
    550493       IF ( .NOT. aero_resist_kray )  THEN
    551           IF ( most_method == 'circular' )  THEN
    552              DO  l = 0, 1
    553                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    554                    surf => surf_lsm_v(l)
    555 !
    556 !--                Compute scaling parameters
    557                    CALL calc_scaling_parameters
    558 !
    559 !--                Compute surface-parallel velocity
    560                    CALL calc_uvw_abs_v_ugrid
    561 !
    562 !--                Compute Obukhov length
    563                    IF ( .NOT. neutral )  CALL calc_ol
    564 !
    565 !--                Compute respective friction velocity on staggered grid
    566                    CALL calc_us
    567 !
    568 !--                Compute respective surface fluxes for momentum and TKE
    569                    CALL calc_surface_fluxes
    570                 ENDIF
    571              ENDDO
    572           ELSE
    573              DO  l = 0, 1
    574                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    575                    surf => surf_lsm_v(l)
    576 !
    577 !--                Compute surface-parallel velocity
    578                    CALL calc_uvw_abs_v_ugrid
    579 !
    580 !--                Compute Obukhov length
    581                    IF ( .NOT. neutral )  CALL calc_ol
    582 !
    583 !--                Compute respective friction velocity on staggered grid
    584                    CALL calc_us
    585 !
    586 !--                Compute scaling parameters
    587                    CALL calc_scaling_parameters
    588 !
    589 !--                Compute respective surface fluxes for momentum and TKE
    590                    CALL calc_surface_fluxes
    591                 ENDIF
    592              ENDDO
    593           ENDIF
     494          DO  l = 0, 1
     495             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     496                surf => surf_lsm_v(l)
     497!
     498!--             Compute surface-parallel velocity
     499                CALL calc_uvw_abs_v_ugrid
     500!
     501!--             Compute Obukhov length
     502                IF ( .NOT. neutral )  CALL calc_ol
     503!
     504!--             Compute respective friction velocity on staggered grid
     505                CALL calc_us
     506!
     507!--             Compute scaling parameters
     508                CALL calc_scaling_parameters
     509!
     510!--             Compute respective surface fluxes for momentum and TKE
     511                CALL calc_surface_fluxes
     512             ENDIF
     513          ENDDO
    594514!
    595515!--    No ts is required, so scaling parameters and Obukhov length do not need
     
    652572!--    flux in land-surface model in case of aero_resist_kray is not true.
    653573       IF ( .NOT. aero_resist_kray )  THEN
    654           IF ( most_method == 'circular' )  THEN
    655              DO  l = 2, 3
    656                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    657                    surf => surf_lsm_v(l)
    658 !
    659 !--                Compute scaling parameters
    660                    CALL calc_scaling_parameters
    661 !
    662 !--                Compute surface-parallel velocity
    663                    CALL calc_uvw_abs_v_vgrid
    664 !
    665 !--                Compute Obukhov length
    666                    IF ( .NOT. neutral )  CALL calc_ol
    667 !
    668 !--                Compute respective friction velocity on staggered grid
    669                    CALL calc_us
    670 !
    671 !--                Compute respective surface fluxes for momentum and TKE
    672                    CALL calc_surface_fluxes
    673                 ENDIF
    674              ENDDO
    675           ELSE
    676              DO  l = 2, 3
    677                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    678                    surf => surf_lsm_v(l)
    679 !
    680 !--                Compute surface-parallel velocity
    681                    CALL calc_uvw_abs_v_vgrid
    682 !
    683 !--                Compute Obukhov length
    684                    IF ( .NOT. neutral )  CALL calc_ol
    685 !
    686 !--                Compute respective friction velocity on staggered grid
    687                    CALL calc_us
    688 !
    689 !--                Compute scaling parameters
    690                    CALL calc_scaling_parameters
    691 !
    692 !--                Compute respective surface fluxes for momentum and TKE
    693                    CALL calc_surface_fluxes
    694                 ENDIF
    695              ENDDO
    696           ENDIF
     574          DO  l = 2, 3
     575             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     576                surf => surf_lsm_v(l)
     577!
     578!--             Compute surface-parallel velocity
     579                CALL calc_uvw_abs_v_vgrid
     580!
     581!--             Compute Obukhov length
     582                IF ( .NOT. neutral )  CALL calc_ol
     583!
     584!--             Compute respective friction velocity on staggered grid
     585                CALL calc_us
     586!
     587!--             Compute scaling parameters
     588                CALL calc_scaling_parameters
     589!
     590!--             Compute respective surface fluxes for momentum and TKE
     591                CALL calc_surface_fluxes
     592             ENDIF
     593          ENDDO
    697594       ELSE
    698595          DO  l = 2, 3
     
    841738          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
    842739          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
    843        ENDIF
    844 
    845        IF ( most_method == 'lookup' )  THEN
    846 
    847 !
    848 !--       Check for roughness heterogeneity. In that case terminate run and
    849 !--       inform user. Check for both, natural and non-natural walls.
    850           IF ( surf_def_h(0)%ns >= 1 )  THEN
    851              IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR. &
    852                   MINVAL( surf_def_h(0)%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
    853                 terminate_run_l = .TRUE.
    854              ENDIF
    855           ENDIF
    856           IF ( surf_lsm_h%ns >= 1 )  THEN
    857              IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
    858                   MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
    859                 terminate_run_l = .TRUE.
    860              ENDIF
    861           ENDIF
    862           IF ( surf_usm_h%ns >= 1 )  THEN
    863              IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h )  .OR.       &
    864                   MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_usm_h%z0  ) )  THEN
    865                 terminate_run_l = .TRUE.
    866              ENDIF
    867           ENDIF
    868 !
    869 !--       Check roughness homogeneity between differt surface types.
    870           IF ( surf_lsm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
    871              IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
    872                   MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
    873                 terminate_run_l = .TRUE.
    874              ENDIF
    875           ENDIF
    876           IF ( surf_usm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
    877              IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
    878                   MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
    879                 terminate_run_l = .TRUE.
    880              ENDIF
    881           ENDIF
    882           IF ( surf_usm_h%ns >= 1  .AND.  surf_lsm_h%ns >= 1 )  THEN
    883              IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
    884                   MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
    885                 terminate_run_l = .TRUE.
    886              ENDIF
    887           ENDIF
    888 
    889 
    890 #if defined( __parallel )
    891 !
    892 !--       Make a logical OR for all processes. Force termiation of model if result
    893 !--       is TRUE
    894           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    895           CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
    896                               MPI_LOR, comm2d, ierr )
    897 #else
    898           terminate_run = terminate_run_l
    899 #endif
    900 
    901           IF ( terminate_run )  THEN
    902              message_string = 'most_method = "lookup" cannot be used in ' //   &
    903                               'combination with a prescribed roughness '  //   &
    904                               'heterogeneity'
    905              CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )
    906           ENDIF
    907 
    908           ALLOCATE(  zeta_tmp(0:num_steps-1) )
    909 
    910 !
    911 !--       Use the lowest possible value for z_mo
    912           k    = nzb
    913           z_mo = zu(k+1) - zw(k)
    914 
    915 !
    916 !--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
    917 !--       available steps (num_steps). The calculation is done with negative
    918 !--       values of zeta in order to simplify the stretching in the free
    919 !--       convection limit for the remaining 10% of steps.
    920           zeta_tmp(0) = - zeta_max
    921           num_steps_n = ( num_steps * 9 / 10 ) - 1
    922           zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
    923 
    924           DO li = 1, num_steps_n
    925              zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
    926           ENDDO
    927 
    928 !
    929 !--       Calculate stretching factor for the free convection range
    930           DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
    931              regr_old = regr
    932              regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
    933                     )**( 10.0_wp / REAL(num_steps) )
    934           ENDDO
    935 
    936 !
    937 !--       Calculate z/L range from zeta_min to zeta_stretch
    938           DO li = num_steps_n+1, num_steps-1
    939              zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
    940              zeta_step = zeta_step * regr
    941           ENDDO
    942 
    943 !
    944 !--       Save roughness lengths to temporary variables
    945           IF ( surf_def_h(0)%ns >= 1  )  THEN
    946              z0h_min = surf_def_h(0)%z0h(1)
    947              z0_min  = surf_def_h(0)%z0(1)
    948           ELSEIF ( surf_lsm_h%ns >= 1 )  THEN
    949              z0h_min = surf_lsm_h%z0h(1)
    950              z0_min  = surf_lsm_h%z0(1)
    951           ELSEIF ( surf_usm_h%ns >= 1 )  THEN
    952              z0h_min = surf_usm_h%z0h(1)
    953              z0_min  = surf_usm_h%z0(1)
    954           ENDIF         
    955 !
    956 !--       Calculate lookup table for the Richardson number versus Obukhov length
    957 !--       The Richardson number (rib) is defined depending on the choice of
    958 !--       boundary conditions for temperature
    959           IF ( ibc_pt_b == 1 )  THEN
    960              DO li = 0, num_steps-1
    961                 ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
    962                 rib_tab(li) = z_mo / ol_tab(li)  / ( LOG( z_mo / z0_min )       &
    963                                                 - psi_m( z_mo / ol_tab(li) )    &
    964                                                 + psi_m( z0_min / ol_tab(li) )  &
    965                                                   )**3
    966              ENDDO 
    967           ELSE
    968              DO li = 0, num_steps-1
    969                 ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
    970                 rib_tab(li) = z_mo / ol_tab(li)  * ( LOG( z_mo / z0h_min )     &
    971                                               - psi_h( z_mo / ol_tab(li) )     &
    972                                               + psi_h( z0h_min / ol_tab(li) )  &
    973                                             )                                  &
    974                                           / ( LOG( z_mo / z0_min )             &
    975                                               - psi_m( z_mo / ol_tab(li) )     &
    976                                               + psi_m( z0_min / ol_tab(li) )   &
    977                                             )**2
    978              ENDDO
    979           ENDIF
    980 
    981 !
    982 !--       Determine minimum values of rib in the lookup table. Set upper limit
    983 !--       to critical Richardson number (0.25)
    984           rib_min  = MINVAL(rib_tab)
    985           rib_max  = 0.25 !MAXVAL(rib_tab)
    986 
    987           DEALLOCATE( zeta_tmp )
    988740       ENDIF
    989741
     
    12561008                       ol_u      !< Upper bound of L for Newton iteration
    12571009
    1258        IF ( TRIM( most_method ) /= 'circular' )  THEN
    1259 !
    1260 !--       Evaluate bulk Richardson number (calculation depends on
    1261 !--       definition based on setting of boundary conditions
    1262           IF ( ibc_pt_b /= 1 )  THEN
    1263              IF ( humidity )  THEN
    1264                 !$OMP PARALLEL DO PRIVATE( z_mo )
    1265                 DO  m = 1, surf%ns
    1266 
    1267                    z_mo = surf%z_mo(m)
    1268 
    1269                    surf%rib(m) = g * z_mo                                      &
    1270                                  * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
    1271                                  / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
    1272                                  + 1.0E-20_wp )
    1273                 ENDDO
    1274              ELSE
    1275                 !$OMP PARALLEL DO PRIVATE( z_mo )
    1276                 DO  m = 1, surf%ns
    1277 
    1278                    z_mo = surf%z_mo(m)
    1279 
    1280                    surf%rib(m) = g * z_mo                                      &
    1281                                  * ( surf%pt1(m) - surf%pt_surface(m) )        &
    1282                                  / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
    1283                                  + 1.0E-20_wp )
    1284                 ENDDO
    1285              ENDIF
     1010!
     1011!--    Evaluate bulk Richardson number (calculation depends on
     1012!--    definition based on setting of boundary conditions
     1013       IF ( ibc_pt_b /= 1 )  THEN
     1014          IF ( humidity )  THEN
     1015             !$OMP PARALLEL DO PRIVATE( z_mo )
     1016             DO  m = 1, surf%ns
     1017
     1018                z_mo = surf%z_mo(m)
     1019
     1020                surf%rib(m) = g * z_mo                                         &
     1021                              * ( surf%vpt1(m) - surf%vpt_surface(m) )         &
     1022                              / ( surf%uvw_abs(m)**2 * surf%vpt1(m)            &
     1023                              + 1.0E-20_wp )
     1024             ENDDO
    12861025          ELSE
    1287              IF ( humidity )  THEN
    1288                 !$OMP PARALLEL DO PRIVATE( k, z_mo )
    1289                 DO  m = 1, surf%ns
    1290 
    1291                    k   = surf%k(m)
    1292 
    1293                    z_mo = surf%z_mo(m)
    1294 
    1295                    surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
     1026             !$OMP PARALLEL DO PRIVATE( z_mo )
     1027             DO  m = 1, surf%ns
     1028
     1029                z_mo = surf%z_mo(m)
     1030
     1031                surf%rib(m) = g * z_mo                                         &
     1032                              * ( surf%pt1(m) - surf%pt_surface(m) )           &
     1033                              / ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp )
     1034             ENDDO
     1035          ENDIF
     1036       ELSE
     1037          IF ( humidity )  THEN
     1038             !$OMP PARALLEL DO PRIVATE( k, z_mo )
     1039             DO  m = 1, surf%ns
     1040
     1041                k   = surf%k(m)
     1042
     1043                z_mo = surf%z_mo(m)
     1044
     1045                surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp                &
    12961046                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
    12971047                           * surf%pt1(m) * surf%qsws(m) ) *                    &
     
    12991049                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
    13001050                           + 1.0E-20_wp )
    1301                 ENDDO
     1051             ENDDO
     1052          ELSE
     1053             !$OMP PARALLEL DO PRIVATE( k, z_mo )
     1054             !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
     1055             !$ACC PRESENT(surf, drho_air_zw)
     1056             DO  m = 1, surf%ns
     1057
     1058                k   = surf%k(m)
     1059
     1060                z_mo = surf%z_mo(m)
     1061
     1062                surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
     1063                                     drho_air_zw(k-1)            /          &
     1064                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2       &
     1065                        + 1.0E-20_wp )
     1066             ENDDO
     1067          ENDIF
     1068       ENDIF
     1069
     1070!
     1071!--    Calculate the Obukhov length using Newton iteration
     1072       !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
     1073       !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
     1074       !$ACC PRESENT(surf)
     1075       DO  m = 1, surf%ns
     1076
     1077          i   = surf%i(m)           
     1078          j   = surf%j(m)
     1079
     1080          z_mo = surf%z_mo(m)
     1081
     1082!
     1083!--       Store current value in case the Newton iteration fails
     1084          ol_old = surf%ol(m)
     1085
     1086!
     1087!--       Ensure that the bulk Richardson number and the Obukhov
     1088!--       length have the same sign
     1089          IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
     1090               ABS( surf%ol(m) ) == ol_max )  THEN
     1091             IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
     1092             IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
     1093          ENDIF
     1094!
     1095!--       Iteration to find Obukhov length
     1096          iter = 0
     1097          DO
     1098             iter = iter + 1
     1099!
     1100!--          In case of divergence, use the value of the previous time step
     1101             IF ( iter > 1000 )  THEN
     1102                surf%ol(m) = ol_old
     1103                EXIT
     1104             ENDIF
     1105
     1106             ol_m = surf%ol(m)
     1107             ol_l = ol_m - 0.001_wp * ol_m
     1108             ol_u = ol_m + 0.001_wp * ol_m
     1109
     1110
     1111             IF ( ibc_pt_b /= 1 )  THEN
     1112!
     1113!--             Calculate f = Ri - [...]/[...]^2 = 0
     1114                f = surf%rib(m) - ( z_mo / ol_m ) * (                          &
     1115                                              LOG( z_mo / surf%z0h(m) )        &
     1116                                              - psi_h( z_mo / ol_m )           &
     1117                                              + psi_h( surf%z0h(m) /           &
     1118                                                       ol_m )                  & 
     1119                                                     )                         &
     1120                                           / ( LOG( z_mo / surf%z0(m) )        &
     1121                                              - psi_m( z_mo / ol_m )           &
     1122                                              + psi_m( surf%z0(m) /  ol_m )    &
     1123                                                 )**2
     1124
     1125!
     1126!--              Calculate df/dL
     1127                 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /                  &
     1128                                                          surf%z0h(m) )        &
     1129                                         - psi_h( z_mo / ol_u )                &
     1130                                         + psi_h( surf%z0h(m) / ol_u )         &
     1131                                           )                                   &
     1132                                         / ( LOG( z_mo / surf%z0(m) )          &
     1133                                         - psi_m( z_mo / ol_u )                &
     1134                                         + psi_m( surf%z0(m) / ol_u )          &
     1135                                           )**2                                &
     1136                        + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )        &
     1137                                         - psi_h( z_mo / ol_l )                &
     1138                                         + psi_h( surf%z0h(m) / ol_l )         &
     1139                                           )                                   &
     1140                                         / ( LOG( z_mo / surf%z0(m) )          &
     1141                                         - psi_m( z_mo / ol_l )                &
     1142                                         + psi_m( surf%z0(m) / ol_l )          &
     1143                                           )**2                                &
     1144                          ) / ( ol_u - ol_l )
    13021145             ELSE
    1303                 !$OMP PARALLEL DO PRIVATE( k, z_mo )
    1304                 !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
    1305                 !$ACC PRESENT(surf, drho_air_zw)
    1306                 DO  m = 1, surf%ns
    1307 
    1308                    k   = surf%k(m)
    1309 
    1310                    z_mo = surf%z_mo(m)
    1311 
    1312                    surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
    1313                                         drho_air_zw(k-1)            /          &
    1314                         ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
    1315                            + 1.0E-20_wp )
    1316                 ENDDO
     1146!
     1147!--             Calculate f = Ri - 1 /[...]^3 = 0
     1148                f = surf%rib(m) - ( z_mo / ol_m ) /                            &
     1149                                             ( LOG( z_mo / surf%z0(m) )        &
     1150                                         - psi_m( z_mo / ol_m )                &
     1151                                         + psi_m( surf%z0(m) / ol_m )          &
     1152                                             )**3
     1153
     1154!
     1155!--             Calculate df/dL
     1156                f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )      &
     1157                                         - psi_m( z_mo / ol_u )                &
     1158                                         + psi_m( surf%z0(m) / ol_u )          &
     1159                                                  )**3                         &
     1160                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
     1161                                         - psi_m( z_mo / ol_l )                &
     1162                                         + psi_m( surf%z0(m) / ol_l )          &
     1163                                            )**3                               &
     1164                          ) / ( ol_u - ol_l )
    13171165             ENDIF
    1318           ENDIF
    1319 
    1320        ENDIF
    1321 
    1322 
    1323 !
    1324 !--    Calculate the Obukhov length either using a Newton iteration
    1325 !--    method, via a lookup table, or using the old circular way
    1326        IF ( TRIM( most_method ) == 'newton' )  THEN
    1327 
    1328           !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
    1329           !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
    1330           !$ACC PRESENT(surf)
    1331           DO  m = 1, surf%ns
    1332 
    1333              i   = surf%i(m)           
    1334              j   = surf%j(m)
    1335 
    1336              z_mo = surf%z_mo(m)
    1337 
    1338 !
    1339 !--          Store current value in case the Newton iteration fails
    1340              ol_old = surf%ol(m)
     1166!
     1167!--          Calculate new L
     1168             surf%ol(m) = ol_m - f / f_d_ol
    13411169
    13421170!
    13431171!--          Ensure that the bulk Richardson number and the Obukhov
    1344 !--          length have the same sign
    1345              IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
    1346                   ABS( surf%ol(m) ) == ol_max )  THEN
    1347                 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
    1348                 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
     1172!--          length have the same sign and ensure convergence.
     1173             IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
     1174!
     1175!--          If unrealistic value occurs, set L to the maximum
     1176!--          value that is allowed
     1177             IF ( ABS( surf%ol(m) ) > ol_max )  THEN
     1178                surf%ol(m) = ol_max
     1179                EXIT
    13491180             ENDIF
    13501181!
    1351 !--          Iteration to find Obukhov length
    1352              iter = 0
    1353              DO
    1354                 iter = iter + 1
    1355 !
    1356 !--             In case of divergence, use the value of the previous time step
    1357                 IF ( iter > 1000 )  THEN
    1358                    surf%ol(m) = ol_old
    1359                    EXIT
    1360                 ENDIF
    1361 
    1362                 ol_m = surf%ol(m)
    1363                 ol_l = ol_m - 0.001_wp * ol_m
    1364                 ol_u = ol_m + 0.001_wp * ol_m
    1365 
    1366 
    1367                 IF ( ibc_pt_b /= 1 )  THEN
    1368 !
    1369 !--                Calculate f = Ri - [...]/[...]^2 = 0
    1370                    f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
    1371                                                  LOG( z_mo / surf%z0h(m) )     &
    1372                                                  - psi_h( z_mo / ol_m )        &
    1373                                                  + psi_h( surf%z0h(m) /        &
    1374                                                           ol_m )               &
    1375                                                                )               &
    1376                                               / ( LOG( z_mo / surf%z0(m) )     &
    1377                                                  - psi_m( z_mo / ol_m )        &
    1378                                                  + psi_m( surf%z0(m) /         &
    1379                                                           ol_m )               &
    1380                                                  )**2
    1381 
    1382 !
    1383 !--                 Calculate df/dL
    1384                     f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
    1385                                                              surf%z0h(m) )     &
    1386                                             - psi_h( z_mo / ol_u )             &
    1387                                             + psi_h( surf%z0h(m) / ol_u )      &
    1388                                               )                                &
    1389                                             / ( LOG( z_mo / surf%z0(m) )       &
    1390                                             - psi_m( z_mo / ol_u )             &
    1391                                             + psi_m( surf%z0(m) / ol_u )       &
    1392                                               )**2                             &
    1393                            + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
    1394                                             - psi_h( z_mo / ol_l )             &
    1395                                             + psi_h( surf%z0h(m) / ol_l )      &
    1396                                               )                                &
    1397                                             / ( LOG( z_mo / surf%z0(m) )       &
    1398                                             - psi_m( z_mo / ol_l )             &
    1399                                             + psi_m( surf%z0(m) / ol_l )       &
    1400                                               )**2                             &
    1401                              ) / ( ol_u - ol_l )
    1402                 ELSE
    1403 !
    1404 !--                Calculate f = Ri - 1 /[...]^3 = 0
    1405                    f = surf%rib(m) - ( z_mo / ol_m ) /                         &
    1406                                                 ( LOG( z_mo / surf%z0(m) )     &
    1407                                             - psi_m( z_mo / ol_m )             &
    1408                                             + psi_m( surf%z0(m) / ol_m )       &
    1409                                                 )**3
    1410 
    1411 !
    1412 !--                Calculate df/dL
    1413                    f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
    1414                                                             surf%z0(m) )       &
    1415                                             - psi_m( z_mo / ol_u )             &
    1416                                             + psi_m( surf%z0(m) / ol_u )       &
    1417                                                      )**3                      &
    1418                            + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
    1419                                             - psi_m( z_mo / ol_l )             &
    1420                                             + psi_m( surf%z0(m) / ol_l )       &
    1421                                                )**3                            &
    1422                                   ) / ( ol_u - ol_l )
    1423                 ENDIF
    1424 !
    1425 !--             Calculate new L
    1426                 surf%ol(m) = ol_m - f / f_d_ol
    1427 
    1428 !
    1429 !--             Ensure that the bulk Richardson number and the Obukhov
    1430 !--             length have the same sign and ensure convergence.
    1431                 IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
    1432 !
    1433 !--             If unrealistic value occurs, set L to the maximum
    1434 !--             value that is allowed
    1435                 IF ( ABS( surf%ol(m) ) > ol_max )  THEN
    1436                    surf%ol(m) = ol_max
    1437                    EXIT
    1438                 ENDIF
    1439 !
    1440 !--             Check for convergence
    1441                 IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
    1442                      surf%ol(m) ) < 1.0E-4_wp )  THEN
    1443                    EXIT
    1444                 ELSE
    1445                    CYCLE
    1446                 ENDIF
    1447 
    1448              ENDDO
     1182!--          Check for convergence
     1183             IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
     1184                EXIT
     1185             ELSE
     1186                CYCLE
     1187             ENDIF
     1188
    14491189          ENDDO
    1450 
    1451        ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
    1452 
    1453           !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
    1454           DO  m = 1, surf%ns
    1455 
    1456              i   = surf%i(m)           
    1457              j   = surf%j(m)
    1458 !
    1459 !--          If the bulk Richardson number is outside the range of the lookup
    1460 !--          table, set it to the exceeding threshold value
    1461              IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
    1462              IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
    1463 
    1464 !
    1465 !--          Find the correct index bounds for linear interpolation. As the
    1466 !--          Richardson number will not differ very much from time step to
    1467 !--          time step , use the index from the last step and search in the
    1468 !--          correct direction
    1469              li = li_bnd
    1470              IF ( rib_tab(li) - surf%rib(m) > 0.0_wp )  THEN
    1471                 DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp  .AND.  li > 0 )
    1472                    li = li-1
    1473                 ENDDO
    1474              ELSE
    1475                 DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
    1476                            .AND.  li < num_steps-1 )
    1477                    li = li+1
    1478                 ENDDO
    1479              ENDIF
    1480              li_bnd = li
    1481 
    1482 !
    1483 !--          Linear interpolation to find the correct value of z/L
    1484              surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
    1485                          / (  rib_tab(li) - rib_tab(li-1) )                    &
    1486                          * ( surf%rib(m) - rib_tab(li-1) ) )
    1487 
    1488           ENDDO
    1489 
    1490        ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
    1491 
    1492           IF ( .NOT. humidity )  THEN
    1493              !$OMP PARALLEL DO PRIVATE( z_mo )
    1494              DO  m = 1, surf%ns
    1495 
    1496                 z_mo = surf%z_mo(m)
    1497 
    1498                 surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
    1499                                        ( kappa * g *                           &
    1500                                          surf%ts(m) + 1E-30_wp )
    1501 !
    1502 !--             Limit the value range of the Obukhov length.
    1503 !--             This is necessary for very small velocities (u,v --> 1), because
    1504 !--             the absolute value of ol can then become very small, which in
    1505 !--             consequence would result in very large shear stresses and very
    1506 !--             small momentum fluxes (both are generally unrealistic).
    1507                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
    1508                    surf%ol(m) = z_mo / zeta_min
    1509                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
    1510                    surf%ol(m) = z_mo / zeta_max
    1511 
    1512              ENDDO
    1513           ELSEIF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
    1514              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
    1515              DO  m = 1, surf%ns
    1516 
    1517                 i   = surf%i(m)           
    1518                 j   = surf%j(m)
    1519                 k   = surf%k(m)
    1520 
    1521                 z_mo = surf%z_mo(m)
    1522 
    1523 
    1524                 surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
    1525                     ( kappa * g * ( surf%ts(m) +                               &
    1526                          0.61_wp * surf%pt1(m) * surf%us(m)                    &
    1527                        + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
    1528                                    surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
    1529 !
    1530 !--             Limit the value range of the Obukhov length.
    1531 !--             This is necessary for very small velocities (u,v --> 1), because
    1532 !--             the absolute value of ol can then become very small, which in
    1533 !--             consequence would result in very large shear stresses and very
    1534 !--             small momentum fluxes (both are generally unrealistic).
    1535                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
    1536                    surf%ol(m) = z_mo / zeta_min
    1537                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
    1538                    surf%ol(m) = z_mo / zeta_max
    1539 
    1540              ENDDO
    1541           ELSE
    1542 
    1543              !$OMP PARALLEL DO PRIVATE( z_mo )
    1544              DO  m = 1, surf%ns
    1545 
    1546                 z_mo = surf%z_mo(m)
    1547 
    1548                 surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
    1549                     ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
    1550                                     surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
    1551                                     surf%ts(m) ) + 1E-30_wp )
    1552 
    1553 !
    1554 !--             Limit the value range of the Obukhov length.
    1555 !--             This is necessary for very small velocities (u,v --> 1), because
    1556 !--             the absolute value of ol can then become very small, which in
    1557 !--             consequence would result in very large shear stresses and very
    1558 !--             small momentum fluxes (both are generally unrealistic).
    1559                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
    1560                    surf%ol(m) = z_mo / zeta_min
    1561                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
    1562                    surf%ol(m) = z_mo / zeta_max
    1563 
    1564              ENDDO
    1565 
    1566           ENDIF
    1567 
    1568        ENDIF
     1190       ENDDO
    15691191
    15701192    END SUBROUTINE calc_ol
Note: See TracChangeset for help on using the changeset viewer.