Changeset 4757 for palm/trunk/SOURCE


Ignore:
Timestamp:
Oct 26, 2020 10:23:38 AM (4 years ago)
Author:
schwenkel
Message:

implement relative humidity as diagnostic output quantity

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r4671 r4757  
    2424! -----------------
    2525! $Id$
     26! Implement relative humidity as diagnostic output quantity
     27!
     28! 4671 2020-09-09 20:27:58Z pavelkrc
    2629! Implementation of downward facing USM and LSM surfaces
    2730!
     
    110113    USE arrays_3d,                                                                                 &
    111114        ONLY:  ddzu,                                                                               &
     115               exner,                                                                              &
     116               hyp,                                                                                &
    112117               pt,                                                                                 &
    113118               q,                                                                                  &
     
    119124
    120125    USE basic_constants_and_equations_mod,                                                         &
    121         ONLY:  kappa, pi
     126        ONLY:  kappa,                                                                              &
     127               pi,                                                                                 &
     128               magnus,                                                                             &
     129               rd_d_rv
    122130
    123131    USE control_parameters,                                                                        &
    124132        ONLY:  current_timestep_number,                                                            &
    125133               data_output,                                                                        &
     134               data_output_pr,                                                                     &
     135               humidity,                                                                           &
    126136               message_string,                                                                     &
    127137               restart_data_format_output,                                                         &
     
    153163    USE kinds
    154164
     165    USE profil_parameter,                                                                          &
     166        ONLY:  dopr_index
     167
    155168    USE restart_data_mpi_io_mod,                                                                   &
    156169        ONLY:  rd_mpi_io_check_array,                                                              &
    157170               rrd_mpi_io,                                                                         &
    158171               wrd_mpi_io
     172
     173    USE statistics,                                                                                &
     174        ONLY:  hom,                                                                                &
     175               pr_palm,                                                                            &
     176               rmask,                                                                              &
     177               statistic_regions,                                                                  &
     178               sums_l
    159179
    160180    USE surface_mod,                                                                               &
     
    164184               surf_usm_h
    165185
    166 
    167186    IMPLICIT NONE
    168187
     
    179198    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m_av !< averaged horizontal wind speed at 10m
    180199
     200    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  rh     !< relative humidity
     201    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  rh_av  !< avg. relative humidity
    181202    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
    182203    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity
     
    220241           doq_calculate,                                                                          &
    221242           doq_check_data_output,                                                                  &
     243           doq_check_data_output_pr,                                                               &
    222244           doq_define_netcdf_grid,                                                                 &
    223245           doq_init,                                                                               &
     
    225247           doq_output_3d,                                                                          &
    226248           doq_output_mask,                                                                        &
     249           doq_statistics,                                                                         &
    227250           doq_rrd_local,                                                                          &
    228251           doq_wrd_local
     
    241264    END INTERFACE doq_check_data_output
    242265
     266    INTERFACE doq_check_data_output_pr
     267       MODULE PROCEDURE doq_check_data_output_pr
     268    END INTERFACE doq_check_data_output_pr
     269
    243270    INTERFACE doq_define_netcdf_grid
    244271       MODULE PROCEDURE doq_define_netcdf_grid
     
    260287       MODULE PROCEDURE doq_init
    261288    END INTERFACE doq_init
     289
     290    INTERFACE doq_statistics
     291       MODULE PROCEDURE doq_statistics
     292    END INTERFACE doq_statistics
    262293
    263294    INTERFACE doq_prepare
     
    274305    END INTERFACE doq_wrd_local
    275306
    276 
    277307 CONTAINS
    278308
     
    300330       SELECT CASE ( TRIM( variable ) )
    301331
     332          CASE ( 'rh' )
     333             IF ( .NOT. ALLOCATED( rh_av ) )  THEN
     334                ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     335             ENDIF
     336             rh_av = 0.0_wp
     337
    302338          CASE ( 'ti' )
    303339             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
     
    384420
    385421       SELECT CASE ( TRIM( variable ) )
     422
     423          CASE ( 'rh' )
     424             IF ( ALLOCATED( rh_av ) )  THEN
     425                DO  i = nxl, nxr
     426                   DO  j = nys, nyn
     427                      DO  k = nzb, nzt+1
     428                         rh_av(k,j,i) = rh_av(k,j,i) + rh(k,j,i)
     429                      ENDDO
     430                   ENDDO
     431                ENDDO
     432             ENDIF
    386433
    387434          CASE ( 'ti' )
     
    523570       SELECT CASE ( TRIM( variable ) )
    524571
     572          CASE ( 'rh' )
     573             IF ( ALLOCATED( rh_av ) )  THEN
     574                DO  i = nxl, nxr
     575                   DO  j = nys, nyn
     576                      DO  k = nzb, nzt+1
     577                         rh_av(k,j,i) = rh_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     578                      ENDDO
     579                   ENDDO
     580                ENDDO
     581             ENDIF
     582
    525583          CASE ( 'ti' )
    526584             IF ( ALLOCATED( ti_av ) )  THEN
     
    688746    SELECT CASE ( TRIM( var ) )
    689747
     748       CASE ( 'rh' )
     749          unit = '%'
     750
    690751       CASE ( 'ti' )
    691752          unit = '1/s'
     
    740801 END SUBROUTINE doq_check_data_output
    741802
     803
     804!--------------------------------------------------------------------------------------------------!
     805! Description:
     806! ------------
     807!> Set the unit of user defined profile output quantities. For those variables not recognized by the
     808!> user, the parameter unit is set to "illegal", which tells the calling routine that the
     809!> output variable is not defined and leads to a program abort.
     810!--------------------------------------------------------------------------------------------------!
     811 SUBROUTINE doq_check_data_output_pr( variable, var_count, unit, dopr_unit )
     812
     813    CHARACTER (LEN=*) ::  unit      !<
     814    CHARACTER (LEN=*) ::  variable  !<
     815    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
     816
     817    INTEGER(iwp) ::  pr_index       !<
     818    INTEGER(iwp) ::  var_count      !<
     819
     820    SELECT CASE ( TRIM( variable ) )
     821
     822       CASE ( 'rh' )
     823          IF ( .NOT.  humidity )  THEN
     824             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
     825                              ' requires humidity'
     826             CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
     827          ENDIF
     828          pr_index = 130
     829          dopr_index(var_count) = pr_index
     830          dopr_unit     = '%'
     831          unit = dopr_unit
     832          hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
     833
     834       CASE DEFAULT
     835          unit = 'illegal'
     836
     837    END SELECT
     838
     839
     840 END SUBROUTINE doq_check_data_output_pr
     841
     842
    742843!--------------------------------------------------------------------------------------------------!
    743844!
     
    763864!
    764865!--    s grid
    765        CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
     866       CASE ( 'rh', 'rh_xy', 'rh_xz', 'rh_yz',                                                     &
     867              'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
    766868              'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                                     &
    767869              'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' )
     
    854956
    855957    SELECT CASE ( TRIM( variable ) )
     958
     959       CASE ( 'rh_xy', 'rh_xz', 'rh_yz' )
     960           IF ( av == 0 )  THEN
     961              to_be_resorted => rh
     962           ELSE
     963              IF ( .NOT. ALLOCATED( rh_av ) )  THEN
     964                 ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     965                 rh_av = REAL( fill_value, KIND = wp )
     966              ENDIF
     967              to_be_resorted => rh_av
     968           ENDIF
     969           flag_nr = 0
     970
     971           IF ( mode == 'xy' )  grid = 'zu'
    856972
    857973       CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
     
    10951211    SELECT CASE ( TRIM( variable ) )
    10961212
     1213       CASE ( 'rh' )
     1214          IF ( av == 0 )  THEN
     1215             to_be_resorted => rh
     1216          ELSE
     1217             IF ( .NOT. ALLOCATED( rh ) )  THEN
     1218                ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1219                rh = REAL( fill_value, KIND = wp )
     1220             ENDIF
     1221             to_be_resorted => rh
     1222          ENDIF
     1223          flag_nr = 0
     1224
    10971225       CASE ( 'ti' )
    10981226          IF ( av == 0 )  THEN
     
    12801408
    12811409    SELECT CASE ( TRIM( variable ) )
     1410
     1411       CASE ( 'rh' )
     1412          IF ( av == 0 )  THEN
     1413             to_be_resorted => rh
     1414          ELSE
     1415             to_be_resorted => rh_av
     1416          ENDIF
     1417          grid = 's'
     1418          flag_nr = 0
    12821419
    12831420       CASE ( 'ti' )
     
    14501587       SELECT CASE ( TRIM( do_all(ivar) ) )
    14511588!
     1589!--       Allocate array for 'relative humidity'
     1590          CASE ( 'rh' )
     1591             IF ( .NOT. ALLOCATED( rh ) )  THEN
     1592                ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1593                rh = 0.0_wp
     1594             ENDIF
     1595!
    14521596!--       Allocate array for 'turbulence intensity'
    14531597          CASE ( 'ti' )
     
    15561700! Description:
    15571701! ------------
     1702!> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series.
     1703!> This is called for every statistic region sr, but at least for the region "total domain" (sr=0).
     1704!--------------------------------------------------------------------------------------------------!
     1705 SUBROUTINE doq_statistics( mode, sr, tn )
     1706
     1707
     1708    CHARACTER (LEN=*) ::  mode   !<
     1709
     1710!     INTEGER(iwp) ::  i    !<
     1711!     INTEGER(iwp) ::  j    !<
     1712!     INTEGER(iwp) ::  k    !<
     1713    INTEGER(iwp) ::  sr   !<
     1714    INTEGER(iwp) ::  tn   !<
     1715!
     1716!-- Next line is to avoid compiler warning about unused variables. Please remove.
     1717    IF ( sr == 0  .OR.  tn == 0 )  CONTINUE
     1718
     1719    IF ( mode == 'profiles' )  THEN
     1720
     1721    ELSEIF ( mode == 'time_series' )  THEN
     1722
     1723    ENDIF
     1724
     1725 END SUBROUTINE doq_statistics
     1726
     1727
     1728!--------------------------------------------------------------------------------------------------!
     1729! Description:
     1730! ------------
    15581731!> Calculate diagnostic quantities
    15591732!--------------------------------------------------------------------------------------------------!
     
    15671740    INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
    15681741
     1742    REAL(wp) ::  temp           !< temperature
     1743    REAL(wp) ::  e_s            !< saturation vapor pressure
     1744    REAL(wp) ::  q_s            !< saturation mixing ratio
     1745
    15691746    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
    15701747
     
    15831760
    15841761       SELECT CASE ( TRIM( do_all(ivar) ) )
     1762!
     1763!--       rh (relative humidity)
     1764          CASE ( 'rh' )
     1765             DO  i = nxl, nxr
     1766                DO  j = nys, nyn
     1767                   DO  k = nzb+1, nzt
     1768                      IF ( humidity )  THEN
     1769                         temp = exner(k) * pt(k,j,i)
     1770                         e_s = magnus( temp )
     1771                         q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
     1772                         rh(k,j,i) = q(k,j,i) / q_s * 100.0_wp                                     &
     1773                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
     1774                      ENDIF
     1775                   ENDDO
     1776                ENDDO
     1777             ENDDO
    15851778!
    15861779!--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
     
    19962189                      tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    19972190
     2191       CASE ( 'rh_av' )
     2192          IF ( .NOT. ALLOCATED( rh_av ) )  THEN
     2193             ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2194          ENDIF
     2195          IF ( k == 1 )  READ ( 13 )  tmp_3d
     2196          rh_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
     2197                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2198
    19982199       CASE ( 'ti_av' )
    19992200          IF ( .NOT. ALLOCATED( ti_av ) )  THEN
     
    21232324    ENDIF
    21242325
     2326    CALL rd_mpi_io_check_array( 'rh_av' , found = array_found )
     2327    IF ( array_found )  THEN
     2328       IF ( .NOT. ALLOCATED( rh_av ) )  ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2329       CALL rrd_mpi_io( 'rh_av', rh_av )
     2330    ENDIF
     2331
    21252332    CALL rd_mpi_io_check_array( 'u_center_av' , found = array_found )
    21262333    IF ( array_found )  THEN
     
    22082415       ENDIF
    22092416
     2417       IF ( ALLOCATED( rh_av ) )  THEN
     2418          CALL wrd_write_string( 'rh_av' )
     2419          WRITE ( 14 )  rh_av
     2420       ENDIF
     2421
    22102422       IF ( ALLOCATED( ti_av ) )  THEN
    22112423          CALL wrd_write_string( 'ti_av' )
     
    22712483
    22722484       IF ( ALLOCATED( pt_2m_av ) )     CALL wrd_mpi_io( 'pt_2m_av', pt_2m_av )
     2485       IF ( ALLOCATED( rh_av ) )        CALL wrd_mpi_io( 'rh_av', rh_av )
    22732486       IF ( ALLOCATED( ti_av ) )        CALL wrd_mpi_io( 'ti_av', ti_av )
    22742487       IF ( ALLOCATED( u_center_av ) )  CALL wrd_mpi_io( 'u_center_av', u_center_av )
  • palm/trunk/SOURCE/dynamics_mod.f90

    r4731 r4757  
    2424! -----------------
    2525! $Id$
     26! Implement relative humidity as diagnostic output quantity
     27!
     28! 4731 2020-10-07 13:25:11Z schwenkel
    2629! Move exchange_horiz from time_integration to modules
    2730!
     
    8285               u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s,                               &
    8386               v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s,                               &
    84                w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s
     87               w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s,                                       &
     88               zu
    8589
    8690    USE basic_constants_and_equations_mod,                                                         &
     
    518522    REAL(wp)     ::  q_s !< saturation mixing ratio
    519523    REAL(wp)     ::  t_l !< actual temperature
     524    REAL(wp)     ::  rh     = 9999999.9_wp !< relative humidity
     525    REAL(wp)     ::  rh_min = 9999999.9_wp !< max relative humidity
     526    REAL(wp)     ::  height = 9999999.9_wp !< height of supersaturated regions
     527    REAL(wp)     ::  min_height = 9999999.9_wp !< height of supersaturated regions
    520528
    521529!
     
    537545                q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
    538546
    539                 IF ( q(k,j,i) > 1.02_wp * q_s )  realistic_q = .FALSE.
     547                IF ( q(k,j,i) > 1.02_wp * q_s )  THEN
     548                   realistic_q = .FALSE.
     549                   rh = q(k,j,i) / q_s * 100.0_wp
     550                   height = zu(k)
     551                ENDIF
    540552             ENDDO
    541553          ENDDO
     
    546558#if defined( __parallel )
    547559       CALL MPI_ALLREDUCE( MPI_IN_PLACE, realistic_q, 1, MPI_LOGICAL, MPI_LAND, comm2d, ierr)
     560       CALL MPI_ALLREDUCE( rh, rh_min, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
     561       CALL MPI_ALLREDUCE( height, min_height, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
    548562#endif
    549563
    550564       IF ( .NOT. realistic_q  .AND.                                                               &
    551565            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    552           message_string = 'The initial mixing ratio exceeds the saturation mixing ratio.'
     566            WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // &
     567               'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
    553568          CALL message( 'dynamic_init_checks', 'PA0697', 2, 2, 0, 6, 0 )
    554569       ELSEIF ( .NOT. realistic_q  .AND.                                                           &
    555570                TRIM( initializing_actions ) == 'read_restart_data' )  THEN
    556           message_string = 'The mixing ratio exceeds the saturation mixing ratio.'
     571          WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' //   &
     572              'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
    557573          CALL message( 'dynamic_init_checks', 'PA0697', 0, 1, 0, 6, 0 )
    558574       ENDIF
  • palm/trunk/SOURCE/flow_statistics.f90

    r4742 r4757  
    2424! -----------------
    2525! $Id$
     26! Implement relative humidity as diagnostic output quantity
     27!
     28! 4742 2020-10-14 15:11:02Z schwenkel
    2629! Implement snow and graupel (bulk microphysics)
    2730!
     
    105108
    106109    USE arrays_3d,                                                                                 &
    107         ONLY:  ddzu, ddzw, d_exner, e, heatflux_output_conversion, hyp, km, kh,                    &
     110        ONLY:  ddzu, ddzw, d_exner, e, exner, heatflux_output_conversion, hyp, km, kh,             &
    108111               momentumflux_output_conversion,                                                     &
    109112               nc, ni, ng, nr, ns, p, prho, prr, pt, q, qc, qi, qg, ql, qr, qs,                    &
     
    112115
    113116    USE basic_constants_and_equations_mod,                                                         &
    114         ONLY:  g, lv_d_cp
     117        ONLY:  g, lv_d_cp, magnus, rd_d_rv
    115118
    116119    USE bulk_cloud_model_mod,                                                                      &
     
    220223    REAL(wp) ::  v2               !<
    221224    REAL(wp) ::  w2               !<
     225
     226    REAL(wp) ::  temp           !< temperature
     227    REAL(wp) ::  e_s            !< saturation vapor pressure
     228    REAL(wp) ::  q_s            !< saturation mixing ratio
     229    REAL(wp) ::  rh_tmp         !< local value of relative humidity
    222230
    223231    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
     
    13901398                      END IF
    13911399                   END IF
     1400                   temp = exner(k) * pt(k,j,i)
     1401                   e_s = magnus( temp )
     1402                   q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
     1403                   rh_tmp = q(k,j,i) / q_s * 100.0_wp
     1404                   sums_l(k,130,tn) = sums_l(k,130,tn) + rh_tmp * rmask(j,i,sr) * flag
    13921405                ENDIF
    13931406!
     
    17471760                IF ( max_pr_cs > 0 )  THEN
    17481761                     sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) =          &
    1749                                sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,0) + &
    1750                                sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,i)
     1762                               sums_l(:,pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs,0) + &
     1763                               sums_l(:,pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs,i)
    17511764
    17521765                ENDIF
     
    18321845             sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
    18331846             sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
    1834              sums(k,123:129)       = sums(k,123:129) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
     1847             sums(k,123:130)       = sums(k,123:130) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
    18351848          ENDIF
    18361849       ENDDO
     
    19461959       hom(:,1,128,sr) = sums(:,128)   ! ns
    19471960       hom(:,1,129,sr) = sums(:,129)   ! qs
     1961       hom(:,1,130,sr) = sums(:,130)   ! rh
    19481962       hom(:,1,73,sr) = sums(:,73)     ! nr
    19491963       hom(:,1,74,sr) = sums(:,74)     ! qr
  • palm/trunk/SOURCE/module_interface.f90

    r4753 r4757  
    2424! -----------------
    2525! $Id$
     26! Implement doq calls
     27!
     28! 4753 2020-10-21 14:55:41Z raasch
    2629! file re-formatted to follow the PALM coding standard
    2730!
     
    342345        ONLY:  doq_3d_data_averaging,                                                              &
    343346               doq_check_data_output,                                                              &
     347               doq_check_data_output_pr,                                                           &
    344348               doq_define_netcdf_grid,                                                             &
    345349               doq_init,                                                                           &
    346350               doq_output_2d,                                                                      &
    347351               doq_output_3d,                                                                      &
     352               doq_statistics,                                                                     &
    348353               doq_rrd_local,                                                                      &
    349354               doq_wrd_local
     
    870875
    871876    CALL dynamics_check_data_output_pr( variable, var_count, unit, dopr_unit )
     877    CALL doq_check_data_output_pr( variable, var_count, unit, dopr_unit )
    872878
    873879    IF ( unit == 'illegal'  .AND.  bulk_cloud_model )  THEN
     
    17011707
    17021708    CALL dynamics_statistics( mode, sr, tn )
     1709    CALL doq_statistics( mode, sr, tn )
    17031710
    17041711    IF ( gust_module_enabled )  CALL gust_statistics( mode, sr, tn, dots_max )
Note: See TracChangeset for help on using the changeset viewer.