Changeset 4331 for palm/trunk


Ignore:
Timestamp:
Dec 10, 2019 6:25:02 PM (5 years ago)
Author:
suehring
Message:

New diagnostic output for 10-m wind speed; Diagnostic output of 2-m potential temperature moved to diagnostic output

Location:
palm/trunk/SOURCE
Files:
15 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r4309 r4331  
    2525# -----------------
    2626# $Id$
     27# Move diagnostic surface output to diagnostic_output_quantities
     28#
     29# 4309 2019-11-26 18:49:59Z suehring
    2730# Add dependency to parallel random generator for synthetic turbulence generator
    2831#
     
    519522diagnostic_output_quantities_mod.o: \
    520523        mod_kinds.o \
    521         modules.o
     524        modules.o \
     525        surface_layer_fluxes_mod.o
    522526diffusion_s.o: \
    523527        mod_kinds.o \
     
    11281132        basic_constants_and_equations_mod.o \
    11291133        cpulog_mod.o \
     1134        diagnostic_output_quantities_mod.o \
    11301135        land_surface_model_mod.o \
    11311136        mod_kinds.o \
  • palm/trunk/SOURCE/average_3d_data.f90

    r4182 r4331  
    2525! -----------------
    2626! $Id$
     27! Move 2-m potential temperature output to diagnostic_output_quantities
     28!
     29! 4182 2019-08-22 15:20:23Z scharf
    2730! Corrected "Former revisions" section
    2831!
     
    315318             ENDIF
    316319
    317          CASE ( 'theta_2m*' )
    318              IF ( ALLOCATED( pt_2m_av ) ) THEN
    319                 DO  i = nxlg, nxrg
    320                    DO  j = nysg, nyng
    321                       pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
    322                    ENDDO
    323                 ENDDO
    324                 CALL exchange_horiz_2d( pt_2m_av, nbgp )
    325              ENDIF
    326 
    327320          CASE ( 't*' )
    328321             IF ( ALLOCATED( ts_av ) ) THEN
  • palm/trunk/SOURCE/check_parameters.f90

    r4301 r4331  
    2525! -----------------
    2626! 4172 2019-08-20 11:55:33Z oliver.maas
     27! Move 2-m potential temperature output to diagnostic_output_quantities
     28!
     29! 11:55:33Z oliver.maas
    2730! removed message PA0421, concerning old parameter recycling_yshift
    2831!
     
    23852388
    23862389          CASE ( 'ghf*', 'lwp*', 'ol*', 'qsws*', 'r_a*',                       &
    2387                  'shf*', 'ssws*', 't*', 'theta_2m*', 'tsurf*', 'us*',          &
     2390                 'shf*', 'ssws*', 't*', 'tsurf*', 'us*',                       &
    23882391                 'z0*', 'z0h*', 'z0q*' )
    23892392             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     
    24242427                                 'res passive_scalar = .TRUE.'
    24252428                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
    2426              ENDIF
    2427 !
    2428 !--          Activate calculation of 2m temperature if output is requested
    2429              IF ( TRIM( var ) == 'theta_2m*' )  THEN
    2430                 do_output_at_2m = .TRUE.
    2431                 unit = 'K'
    24322429             ENDIF
    24332430
  • palm/trunk/SOURCE/data_output_2d.f90

    r4329 r4331  
    2525! -----------------
    2626! $Id$
     27! Move 2-m potential temperature output to diagnostic_output_quantities
     28!
     29! 4329 2019-12-10 15:46:36Z motisi
    2730! Renamed wall_flags_0 to wall_flags_static_0
    2831!
     
    10011004                ENDIF
    10021005                IF ( mode == 'xy' )  level_z = zu
    1003 
    1004              CASE ( 'theta_2m*_xy' )        ! 2d-array
    1005                 IF ( av == 0 )  THEN
    1006                    DO  m = 1, surf_def_h(0)%ns
    1007                       i = surf_def_h(0)%i(m)
    1008                       j = surf_def_h(0)%j(m)
    1009                       local_pf(i,j,nzb+1) = surf_def_h(0)%pt_2m(m)
    1010                    ENDDO
    1011                    DO  m = 1, surf_lsm_h%ns
    1012                       i = surf_lsm_h%i(m)
    1013                       j = surf_lsm_h%j(m)
    1014                       local_pf(i,j,nzb+1) = surf_lsm_h%pt_2m(m)
    1015                    ENDDO
    1016                    DO  m = 1, surf_usm_h%ns
    1017                       i = surf_usm_h%i(m)
    1018                       j = surf_usm_h%j(m)
    1019                       local_pf(i,j,nzb+1) = surf_usm_h%pt_2m(m)
    1020                    ENDDO
    1021                 ELSE
    1022                    IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
    1023                       ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
    1024                       pt_2m_av = REAL( fill_value, KIND = wp )
    1025                    ENDIF
    1026                    DO  i = nxl, nxr
    1027                       DO  j = nys, nyn
    1028                          local_pf(i,j,nzb+1) = pt_2m_av(j,i)
    1029                       ENDDO
    1030                    ENDDO
    1031                 ENDIF
    1032                 resorted = .TRUE.
    1033                 two_d = .TRUE.
    1034                 level_z(nzb+1) = zu(nzb+1)
    10351006
    10361007             CASE ( 'w_xy', 'w_xz', 'w_yz' )
  • palm/trunk/SOURCE/data_output_mask.f90

    r4329 r4331  
    2525! -----------------
    2626! $Id$
     27! Formatting adjustment
     28!
     29! 4329 2019-12-10 15:46:36Z motisi
    2730! Renamed wall_flags_0 to wall_flags_static_0
    2831!
     
    603606             IF ( .NOT. found )  THEN
    604607                CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf,   &
    605                                    mid)
     608                                      mid)
    606609             ENDIF
    607610!
  • palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90

    r4329 r4331  
    2525! -----------------
    2626! $Id$
     27! - Modularize 2-m potential temperature output
     28! - New output for 10-m wind speed
     29!
     30! 4329 2019-12-10 15:46:36Z motisi
    2731! Renamed wall_flags_0 to wall_flags_static_0
    2832!
     
    7579 
    7680    USE arrays_3d,                                                             &
    77         ONLY:  ddzu, u, v, w
    78 !
    79 !     USE averaging
    80 !
    81 !     USE basic_constants_and_equations_mod,                                     &
    82 !         ONLY:  c_p, lv_d_cp, l_v
    83 !
    84 !     USE bulk_cloud_model_mod,                                                  &
    85 !         ONLY:  bulk_cloud_model
    86 !
     81        ONLY:  ddzu,                                                           &
     82               pt,                                                             &
     83               u,                                                              &
     84               v,                                                              &
     85               w,                                                              &
     86               zu,                                                             &
     87               zw
     88
     89    USE basic_constants_and_equations_mod,                                     &
     90        ONLY:  kappa
     91
    8792    USE control_parameters,                                                    &
    88         ONLY:  current_timestep_number, varnamelength
     93        ONLY:  current_timestep_number,                                        &
     94               data_output,                                                    &
     95               message_string,                                                 &
     96               varnamelength
    8997!
    9098!     USE cpulog,                                                                &
     
    93101   USE grid_variables,                                                         &
    94102        ONLY:  ddx, ddy
    95 !
     103
    96104    USE indices,                                                               &
    97         ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_static_0
    98 !
     105        ONLY:  nbgp,                                                           &
     106               nxl,                                                            &
     107               nxlg,                                                           &
     108               nxr,                                                            &
     109               nxrg,                                                           &
     110               nyn,                                                            &
     111               nyng,                                                           &
     112               nys,                                                            &
     113               nysg,                                                           &
     114               nzb,                                                            &
     115               nzt,                                                            &
     116               wall_flags_static_0
     117
    99118    USE kinds
    100 !
    101 !     USE land_surface_model_mod,                                                &
    102 !         ONLY:  zs
    103 !
    104 !     USE module_interface,                                                      &
    105 !         ONLY:  module_interface_data_output_2d
    106 !
    107 ! #if defined( __netcdf )
    108 !     USE NETCDF
    109 ! #endif
    110 !
    111 !     USE netcdf_interface,                                                      &
    112 !         ONLY:  fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d,       &
    113 !                id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
    114 !                netcdf_data_format, netcdf_handle_error
    115 !
    116 !     USE particle_attributes,                                                   &
    117 !         ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
    118 !                particles, prt_count
    119 !     
    120 !     USE pegrid
    121 !
    122 !     USE surface_mod,                                                           &
    123 !         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
    124 !                surf_lsm_h, surf_usm_h
    125 !
    126 !     USE turbulence_closure_mod,                                                &
    127 !         ONLY:  tcm_data_output_2d
     119
     120    USE surface_mod,                                                           &
     121        ONLY:  surf_def_h,                                                     &
     122               surf_lsm_h,                                                     &
     123               surf_type,                                                      &
     124               surf_usm_h
    128125
    129126
     
    137134    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
    138135
     136    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m     !< 2-m air potential temperature
     137    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m_av  !< averaged 2-m air potential temperature
     138    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m    !< horizontal wind speed at 10m
     139    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m_av !< averaged horizontal wind speed at 10m
     140
    139141    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
    140     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity
    141    
     142    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity   
    142143    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu     !< uu
    143     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av  !< mean of uu
    144    
     144    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av  !< mean of uu   
    145145    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv     !< vv
    146     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av  !< mean of vv
    147    
     146    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av  !< mean of vv   
    148147    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww     !< ww
    149148    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av  !< mean of ww
     
    160159           prepared_diagnostic_output_quantities,                              &
    161160           timestep_number_at_prev_calc,                                       &
     161           pt_2m_av,                                                           &
    162162           ti_av,                                                              &
    163163           uu_av,                                                              &
     164           uv_10m_av,                                                          &
    164165           vv_av,                                                              &
    165166           ww_av
     
    270271             ENDIF
    271272             ww_av = 0.0_wp
     273             
     274          CASE ( 'theta_2m*' )
     275             IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
     276                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
     277             ENDIF
     278             pt_2m_av = 0.0_wp
     279             
     280          CASE ( 'wspeed_10m*' )
     281             IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
     282                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
     283             ENDIF
     284             uv_10m_av = 0.0_wp
    272285
    273286          CASE DEFAULT
     
    323336                ENDDO
    324337             ENDIF
     338             
     339          CASE ( 'theta_2m*' )
     340             IF ( ALLOCATED( pt_2m_av ) ) THEN
     341                DO  i = nxl, nxr
     342                   DO  j = nys, nyn
     343                      pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i)
     344                   ENDDO
     345                ENDDO
     346             ENDIF
     347
     348          CASE ( 'wspeed_10m*' )
     349             IF ( ALLOCATED( uv_10m_av ) ) THEN
     350                DO  i = nxl, nxr
     351                   DO  j = nys, nyn
     352                      uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
     353                   ENDDO
     354                ENDDO
     355             ENDIF
    325356
    326357          CASE DEFAULT
     
    377408             ENDIF
    378409
     410         CASE ( 'theta_2m*' )
     411            IF ( ALLOCATED( pt_2m_av ) ) THEN
     412               DO  i = nxlg, nxrg
     413                  DO  j = nysg, nyng
     414                     pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
     415                  ENDDO
     416               ENDDO
     417               CALL exchange_horiz_2d( pt_2m_av, nbgp )
     418            ENDIF
     419
     420         CASE ( 'wspeed_10m*' )
     421            IF ( ALLOCATED( uv_10m_av ) ) THEN
     422               DO  i = nxlg, nxrg
     423                  DO  j = nysg, nyng
     424                     uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp )
     425                  ENDDO
     426               ENDDO
     427               CALL exchange_horiz_2d( uv_10m_av, nbgp )
     428            ENDIF
     429
    379430       END SELECT
    380431
     
    389440!> Check data output for diagnostic output
    390441!------------------------------------------------------------------------------!
    391  SUBROUTINE doq_check_data_output( var, unit )
     442 SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
    392443
    393444    IMPLICIT NONE
     
    395446    CHARACTER (LEN=*) ::  unit  !<
    396447    CHARACTER (LEN=*) ::  var   !<
     448
     449    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  i     !< Current element of data_output
     450    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  ilen  !< Length of current entry in data_output
     451    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
    397452
    398453    SELECT CASE ( TRIM( var ) )
     
    409464       CASE ( 'ww' )
    410465          unit = 'm2/s2'
    411              
     466!
     467!--    Treat horizotal cross-section output quanatities
     468       CASE ( 'theta_2m*', 'wspeed_10m*' )
     469!
     470!--       Check if output quantity is _xy only.
     471          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     472             message_string = 'illegal value for data_output: "' //            &
     473                              TRIM( var ) // '" & only 2d-horizontal ' //      &
     474                              'cross sections are allowed for this value'
     475             CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
     476          ENDIF
     477
     478          IF ( TRIM( var ) == 'theta_2m*'   )  unit = 'K'
     479          IF ( TRIM( var ) == 'wspeed_10m*' )  unit = 'm/s'
     480
    412481       CASE DEFAULT
    413482          unit = 'illegal'
     
    443512          grid_x = 'x'
    444513          grid_y = 'y'
    445           grid_z = 'zu'   
     514          grid_z = 'zu'
     515!
     516!--    s grid surface variables
     517       CASE ( 'theta_2m*_xy', 'wspeed_10m*' )
     518
     519          grid_x = 'x'
     520          grid_y = 'y'
     521          grid_z = 'zu'
    446522!
    447523!--    u grid
     
    572648         
    573649          IF ( mode == 'xy' )  grid = 'zw'
     650         
     651       CASE ( 'theta_2m*_xy' )        ! 2d-array
     652          IF ( av == 0 )  THEN
     653             DO  i = nxl, nxr
     654                DO  j = nys, nyn
     655                   local_pf(i,j,nzb+1) = pt_2m(j,i)
     656                ENDDO
     657             ENDDO
     658          ELSE
     659             IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
     660                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
     661                pt_2m_av = REAL( fill_value, KIND = wp )
     662             ENDIF
     663             DO  i = nxl, nxr
     664                DO  j = nys, nyn
     665                   local_pf(i,j,nzb+1) = pt_2m_av(j,i)
     666                ENDDO
     667             ENDDO
     668          ENDIF
     669          resorted = .TRUE.
     670          two_d    = .TRUE.
     671          grid     = 'zu1'
     672         
     673       CASE ( 'wspeed_10m*_xy' )        ! 2d-array
     674          IF ( av == 0 )  THEN
     675             DO  i = nxl, nxr
     676                DO  j = nys, nyn
     677                   local_pf(i,j,nzb+1) = uv_10m(j,i)
     678                ENDDO
     679             ENDDO
     680          ELSE
     681             IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
     682                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
     683                uv_10m_av = REAL( fill_value, KIND = wp )
     684             ENDIF
     685             DO  i = nxl, nxr
     686                DO  j = nys, nyn
     687                   local_pf(i,j,nzb+1) = uv_10m_av(j,i)
     688                ENDDO
     689             ENDDO
     690          ENDIF
     691          resorted = .TRUE.
     692          two_d    = .TRUE.
     693          grid     = 'zu1'
    574694
    575695       CASE DEFAULT
     
    776896          flag_nr = 3
    777897
    778 
    779898       CASE DEFAULT
    780899          found = .FALSE.
     
    864983                uu = 0.0_wp
    865984             ENDIF
    866 
    867985!
    868986!--       Allocate array for vv
     
    872990                vv = 0.0_wp
    873991             ENDIF
    874 
    875992!
    876993!--       Allocate array for ww
     
    879996                ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
    880997                ww = 0.0_wp
     998             ENDIF
     999!
     1000!--       Allocate array for 2-m potential temperature
     1001          CASE ( 'theta_2m*' )
     1002             IF ( .NOT. ALLOCATED( pt_2m ) )  THEN
     1003                ALLOCATE( pt_2m(nys:nyn,nxl:nxr) )
     1004                pt_2m = 0.0_wp
     1005             ENDIF
     1006!
     1007!--       Allocate array for 10-m wind speed
     1008          CASE ( 'wspeed_10m*' )
     1009             IF ( .NOT. ALLOCATED( uv_10m ) )  THEN
     1010                ALLOCATE( uv_10m(nys:nyn,nxl:nxr) )
     1011                uv_10m = 0.0_wp
    8811012             ENDIF
    8821013
     
    9001031    IMPLICIT NONE
    9011032
    902     INTEGER(iwp) ::  i, j, k    !< grid loop index in x-, y-, z-direction
     1033    INTEGER(iwp) ::  i          !< grid index x-dimension
     1034    INTEGER(iwp) ::  j          !< grid index y-dimension
     1035    INTEGER(iwp) ::  k          !< grid index z-dimension
    9031036    INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
     1037   
     1038    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
    9041039
    9051040
     
    9731108                ENDDO
    9741109             ENDDO
     1110!
     1111!--       2-m potential temperature
     1112          CASE ( 'theta_2m*' )
     1113!
     1114!--          2-m potential temperature is caluclated from surface arrays. In
     1115!--          case the 2m level is below the first grid point, MOST is applied,
     1116!--          else, linear interpolation between two vertical grid levels is
     1117!--          applied. To access all surfaces, iterate over all horizontally-
     1118!--          upward facing surface types.
     1119             surf => surf_def_h(0)
     1120             CALL calc_pt_2m
     1121             surf => surf_lsm_h
     1122             CALL calc_pt_2m
     1123             surf => surf_usm_h
     1124             CALL calc_pt_2m
     1125!
     1126!--       10-m wind speed
     1127          CASE ( 'wspeed_10m*' )
     1128!
     1129!--          10-m wind speed is caluclated from surface arrays. In
     1130!--          case the 10m level is below the first grid point, MOST is applied,
     1131!--          else, linear interpolation between two vertical grid levels is
     1132!--          applied. To access all surfaces, iterate over all horizontally-
     1133!--          upward facing surface types.
     1134             surf => surf_def_h(0)
     1135             CALL calc_wind_10m
     1136             surf => surf_lsm_h
     1137             CALL calc_wind_10m
     1138             surf => surf_usm_h
     1139             CALL calc_wind_10m
    9751140
    9761141       END SELECT
     
    9811146!     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
    9821147
    983  END SUBROUTINE doq_calculate
    984 
    985 
     1148!
     1149!-- The following block contains subroutines to calculate diagnostic
     1150!-- surface quantities.
     1151    CONTAINS
    9861152!------------------------------------------------------------------------------!
    9871153! Description:
    9881154! ------------
    989 !> ...
     1155!> Calculation of 2-m potential temperature.
     1156!------------------------------------------------------------------------------!
     1157       SUBROUTINE calc_pt_2m
     1158
     1159          USE surface_layer_fluxes_mod,                                        &
     1160              ONLY:  psi_h
     1161
     1162          IMPLICIT NONE
     1163
     1164          INTEGER(iwp) ::  kk     !< running index along the z-dimension
     1165          INTEGER(iwp) ::  m      !< running index for surface elements
     1166 
     1167          DO  m = 1, surf%ns
     1168
     1169             i = surf%i(m)
     1170             j = surf%j(m)
     1171             k = surf%k(m)
     1172!
     1173!--          If 2-m level is below the first grid level, MOST is
     1174!--          used for calculation of 2-m temperature.
     1175             IF ( surf%z_mo(m) > 2.0_wp )  THEN
     1176                pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa           &
     1177                                * ( LOG( 2.0_wp /  surf%z0h(m) )               &
     1178                                  - psi_h( 2.0_wp / surf%ol(m) )               &
     1179                                  + psi_h( surf%z0h(m) / surf%ol(m) ) )
     1180!
     1181!--          If 2-m level is above the first grid level, 2-m temperature
     1182!--          is linearly interpolated between the two nearest vertical grid
     1183!--          levels. Note, since 2-m temperature is only computed for
     1184!--          horizontal upward-facing surfaces, only a vertical
     1185!--          interpolation is necessary.
     1186             ELSE
     1187!
     1188!--             zw(k-1) defines the height of the surface.
     1189                kk = k
     1190                DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp  .AND.  kk <= nzt )
     1191                   kk = kk + 1
     1192                ENDDO
     1193!
     1194!--             kk defines the index of the first grid level >= 2m.
     1195                pt_2m(j,i) = pt(kk-1,j,i) +                                    &
     1196                              ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *            &
     1197                              ( pt(kk,j,i)       - pt(kk-1,j,i) ) /            &
     1198                              ( zu(kk)           - zu(kk-1)     )
     1199             ENDIF
     1200
     1201          ENDDO
     1202
     1203       END SUBROUTINE calc_pt_2m
     1204       
     1205!------------------------------------------------------------------------------!
     1206! Description:
     1207! ------------
     1208!> Calculation of 10-m wind speed.
     1209!------------------------------------------------------------------------------!
     1210       SUBROUTINE calc_wind_10m
     1211
     1212          USE surface_layer_fluxes_mod,                                        &
     1213              ONLY:  psi_m
     1214
     1215          IMPLICIT NONE
     1216
     1217          INTEGER(iwp) ::  kk     !< running index along the z-dimension
     1218          INTEGER(iwp) ::  m      !< running index for surface elements
     1219
     1220          REAL(wp) ::  uv_l !< wind speed at lower grid point
     1221          REAL(wp) ::  uv_u !< wind speed at upper grid point
     1222 
     1223          DO  m = 1, surf%ns
     1224
     1225             i = surf%i(m)
     1226             j = surf%j(m)
     1227             k = surf%k(m)
     1228!
     1229!--          If 10-m level is below the first grid level, MOST is
     1230!--          used for calculation of 10-m temperature.
     1231             IF ( surf%z_mo(m) > 10.0_wp )  THEN
     1232                uv_10m(j,i) = surf%us(m) / kappa                               &
     1233                          * ( LOG( 10.0_wp /  surf%z0(m) )                     &
     1234                              - psi_m( 10.0_wp    / surf%ol(m) )               &
     1235                              + psi_m( surf%z0(m) / surf%ol(m) ) )
     1236!
     1237!--          If 10-m level is above the first grid level, 10-m wind speed
     1238!--          is linearly interpolated between the two nearest vertical grid
     1239!--          levels. Note, since 10-m temperature is only computed for
     1240!--          horizontal upward-facing surfaces, only a vertical
     1241!--          interpolation is necessary.
     1242             ELSE
     1243!
     1244!--             zw(k-1) defines the height of the surface.
     1245                kk = k
     1246                DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp  .AND.  kk <= nzt )
     1247                   kk = kk + 1
     1248                ENDDO
     1249!
     1250!--             kk defines the index of the first grid level >= 10m.
     1251                uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2   &
     1252                           + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
     1253
     1254                uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2   &
     1255                           + ( 0.5_wp * ( v(kk,j,i)   + v(kk,j+1,i)   ) )**2 )
     1256
     1257                uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *        &
     1258                                     ( uv_u              - uv_l     ) /        &
     1259                                     ( zu(kk)            - zu(kk-1) )
     1260
     1261             ENDIF
     1262
     1263          ENDDO
     1264
     1265       END SUBROUTINE calc_wind_10m
     1266
     1267 END SUBROUTINE doq_calculate
     1268
     1269
     1270!------------------------------------------------------------------------------!
     1271! Description:
     1272! ------------
     1273!> Preparation of the diagnostic output, counting of the module-specific
     1274!> output quantities and gathering of the output names.
    9901275!------------------------------------------------------------------------------!
    9911276 SUBROUTINE doq_prepare
     
    11731458    IMPLICIT NONE
    11741459
     1460    IF ( ALLOCATED( pt_2m_av ) )  THEN
     1461       CALL wrd_write_string( 'pt_2m_av' )
     1462       WRITE ( 14 )  pt_2m_av
     1463    ENDIF
     1464
    11751465    IF ( ALLOCATED( ti_av ) )  THEN
    1176        CALL wrd_write_string( 'ti_av' ) 
     1466       CALL wrd_write_string( 'ti_av' )
    11771467       WRITE ( 14 )  ti_av
    11781468    ENDIF
    1179    
     1469
    11801470    IF ( ALLOCATED( uu_av ) )  THEN
    1181        CALL wrd_write_string( 'uu_av' ) 
     1471       CALL wrd_write_string( 'uu_av' )
    11821472       WRITE ( 14 )  uu_av
    11831473    ENDIF
    1184        
     1474
     1475    IF ( ALLOCATED( uv_10m_av ) )  THEN
     1476       CALL wrd_write_string( 'uv_10m_av' )
     1477       WRITE ( 14 )  uv_10m_av
     1478    ENDIF
     1479
    11851480    IF ( ALLOCATED( vv_av ) )  THEN
    1186        CALL wrd_write_string( 'vv_av' ) 
     1481       CALL wrd_write_string( 'vv_av' )
    11871482       WRITE ( 14 )  vv_av
    11881483    ENDIF
    1189        
     1484
    11901485    IF ( ALLOCATED( ww_av ) )  THEN
    1191        CALL wrd_write_string( 'ww_av' ) 
     1486       CALL wrd_write_string( 'ww_av' )
    11921487       WRITE ( 14 )  ww_av
    11931488    ENDIF
  • palm/trunk/SOURCE/module_interface.f90

    r4281 r4331  
    2525! -----------------
    2626! $Id$
     27! Change interface for doq_check_data_output, in order to perform further
     28! output checks.
     29!
     30! 4281 2019-10-29 15:15:39Z schwenkel
    2731! Added dynamics boundary conditions
    2832!
     
    872876
    873877    IF ( unit == 'illegal' )  THEN
    874        CALL doq_check_data_output( variable, unit )
     878       CALL doq_check_data_output( variable, unit, i, ilen, k )
    875879    ENDIF
    876880
  • palm/trunk/SOURCE/modules.f90

    r4329 r4331  
    2525! -----------------
    2626! $Id$
     27! - do_output_at_2m, pt_2m_av
     28!
     29! 4329 2019-12-10 15:46:36Z motisi
    2730! Renamed wall_flags_0 to wall_flags_static_0
    2831!
     
    413416    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av                 !< avg. liquid water path
    414417    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ol_av                  !< avg. Obukhov length
    415     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m_av               !< avg. 2m- air potential temperature
    416418    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qsws_av                !< avg. surface moisture flux
    417419    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  r_a_av                 !< avg. resistance
     
    707709    LOGICAL ::  do2d_at_begin = .FALSE.                          !< namelist parameter
    708710    LOGICAL ::  do3d_at_begin = .FALSE.                          !< namelist parameter
    709     LOGICAL ::  do_output_at_2m = .FALSE.                        !< flag for activating calculation of potential temperature at z = 2 m
    710711    LOGICAL ::  do_sum = .FALSE.                                 !< contribute to time average of profile data?
    711712    LOGICAL ::  dp_external = .FALSE.                            !< namelist parameter
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r4329 r4331  
    2727! -----------------
    2828! $Id$
     29! Typo corrected
     30!
     31! 4329 2019-12-10 15:46:36Z motisi
    2932! Renamed wall_flags_0 to wall_flags_static_0
    3033!
     
    11351138                                        'defined on top of an artificially '// &
    11361139                                        'created building grid point ' //      &
    1137                                         '(emerged from the filering) - ' //    &
     1140                                        '(emerged from the filtering) - ' //   &
    11381141                                        'LAD profile is omitted at this ' //   &
    11391142                                        'grid point: (i,j) = ', i, j
  • palm/trunk/SOURCE/read_restart_data_mod.f90

    r4301 r4331  
    2525! -----------------
    2626! $Id$
     27! Enable restart data for 2-m potential temperature output
     28!
     29! 4301 2019-11-22 12:09:09Z oliver.maas
    2730! removed recycling_yshift
    2831!
     
    101104
    102105    USE diagnostic_output_quantities_mod,                                      &
    103         ONLY:  ti_av, uu_av, vv_av, ww_av
     106        ONLY:  pt_2m_av,                                                       &
     107               ti_av,                                                          &
     108               uu_av,                                                          &
     109               uv_10m_av,                                                      &
     110               vv_av,                                                          &
     111               ww_av
    104112
    105113    USE grid_variables,                                                        &
     
    15251533                   uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
    15261534
     1535                CASE ( 'uv_10m_av' )
     1536                   IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
     1537                      ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
     1538                   ENDIF
     1539                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     1540                   uv_10m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =                           &
     1541                      tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1542
    15271543                CASE ( 'u_m_l' )
    15281544                   IF ( k == 1 )  THEN
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r4182 r4331  
    2525! -----------------
    2626! $Id$
     27! Move 2-m potential temperature output to diagnostic_output_quantities
     28!
     29! 4182 2019-08-22 15:20:23Z scharf
    2730! Corrected "Former revisions" section
    2831!
     
    6467
    6568    USE arrays_3d,                                                             &
    66         ONLY:  dzw, d_exner, e, heatflux_output_conversion, p,    &
    67                pt, q, ql, ql_c, ql_v, s, u, v, vpt, w,                 &
     69        ONLY:  dzw, d_exner, e, heatflux_output_conversion, p,                 &
     70               pt, q, ql, ql_c, ql_v, s, u, v, vpt, w,                         &
    6871               waterflux_output_conversion
    6972
    7073    USE averaging,                                                             &
    7174        ONLY:  e_av, ghf_av, lpt_av, lwp_av, ol_av, p_av, pc_av, pr_av, pt_av, &
    72                pt_2m_av, q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av,     &
     75               q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av,               &
    7376               qv_av, r_a_av, s_av, shf_av, ssws_av, ts_av, tsurf_av, u_av,    &
    7477               us_av, v_av, vpt_av, w_av, z0_av, z0h_av, z0q_av
     
    296299                ENDIF
    297300                vpt_av = 0.0_wp
    298 
    299              CASE ( 'theta_2m*' )
    300                 IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
    301                    ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
    302                 ENDIF
    303                 pt_2m_av = 0.0_wp
    304301
    305302             CASE ( 'w' )
     
    730727             ENDIF
    731728
    732           CASE ( 'theta_2m*' )
    733              IF ( ALLOCATED( pt_2m_av ) ) THEN   
    734                 DO  i = nxl, nxr
    735                    DO  j = nys, nyn
    736                       match_def = surf_def_h(0)%start_index(j,i) <=            &
    737                                   surf_def_h(0)%end_index(j,i)
    738                       match_lsm = surf_lsm_h%start_index(j,i) <=               &
    739                                   surf_lsm_h%end_index(j,i)
    740                       match_usm = surf_usm_h%start_index(j,i) <=               &
    741                                   surf_usm_h%end_index(j,i)
    742 
    743                       IF ( match_def )  THEN
    744                          m = surf_def_h(0)%end_index(j,i)
    745                          pt_2m_av(j,i) = pt_2m_av(j,i) +                       &
    746                                          surf_def_h(0)%pt_2m(m)
    747                       ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
    748                          m = surf_lsm_h%end_index(j,i)
    749                          pt_2m_av(j,i) = pt_2m_av(j,i) +                       &
    750                                          surf_lsm_h%pt_2m(m)
    751                       ELSEIF ( match_usm )  THEN
    752                          m = surf_usm_h%end_index(j,i)
    753                          pt_2m_av(j,i) = pt_2m_av(j,i) +                       &
    754                                          surf_usm_h%pt_2m(m)
    755                       ENDIF
    756                    ENDDO
    757                 ENDDO
    758              ENDIF
    759              
    760              
    761729          CASE ( 't*' )
    762730             IF ( ALLOCATED( ts_av ) ) THEN
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r4298 r4331  
    2626! -----------------
    2727! $Id$
     28! Calculation of diagnostic-only 2-m potential temperature moved to
     29! diagnostic_output_quantities.
     30!
     31! 4298 2019-11-21 15:59:16Z suehring
    2832! Calculation of 2-m temperature adjusted to the case the 2-m level is above
    2933! the first grid point.
     
    108112               constant_waterflux, coupling_mode,                              &
    109113               debug_output_timestep,                                          &
    110                do_output_at_2m, humidity,                                      &
     114               humidity,                                                       &
    111115               ibc_e_b, ibc_pt_b, indoor_model,                                &
    112116               land_surface, large_scale_forcing, lsf_surf, message_string,    &
     
    162166    PRIVATE
    163167
    164     PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
     168    PUBLIC init_surface_layer_fluxes,                                          &
     169           phi_m,                                                              &
     170           psi_h,                                                              &
     171           psi_m,                                                              &
     172           surface_layer_fluxes
    165173
    166174    INTERFACE init_surface_layer_fluxes
     
    171179       MODULE PROCEDURE phi_m
    172180    END INTERFACE phi_m
     181
     182    INTERFACE psi_h
     183       MODULE PROCEDURE psi_h
     184    END INTERFACE psi_h
     185
     186    INTERFACE psi_m
     187       MODULE PROCEDURE psi_m
     188    END INTERFACE psi_m
    173189
    174190    INTERFACE surface_layer_fluxes
     
    259275          CALL calc_scaling_parameters
    260276          CALL calc_surface_fluxes
    261           IF ( do_output_at_2m )  THEN
    262              CALL calc_pt_near_surface ( '2m' )
    263           ENDIF
    264277       ENDIF
    265278!
     
    272285          CALL calc_scaling_parameters
    273286          CALL calc_surface_fluxes
    274           IF ( do_output_at_2m )  THEN
    275              CALL calc_pt_near_surface ( '2m' )
    276           ENDIF
    277287       ENDIF
    278288!
     
    285295          CALL calc_scaling_parameters
    286296          CALL calc_surface_fluxes
    287           IF ( do_output_at_2m )  THEN
    288              CALL calc_pt_near_surface ( '2m' )
    289           ENDIF
    290297!
    291298!--       Calculate 10cm temperature, required in indoor model
     
    18631870       INTEGER(iwp)                    :: j      !< grid index y-dimension
    18641871       INTEGER(iwp)                    :: k      !< grid index z-dimension
    1865        INTEGER(iwp)                    :: kk     !< running index along the z-dimension
    18661872       INTEGER(iwp)                    :: m      !< running index for surface elements
    18671873
     
    18821888                                     - psi_h( 0.1_wp / surf%ol(m) )            &
    18831889                                     + psi_h( surf%z0h(m) / surf%ol(m) ) )
    1884                                    
    1885              ENDDO
    1886 
    1887 !
    1888 !--       2-m temperature. Note, this is only calculated for output reasons at
    1889 !--       horizontal upward-facing surfaces.
    1890           CASE ( '2m' )
    1891      
    1892              DO  m = 1, surf%ns
    1893 
    1894                 i = surf%i(m)
    1895                 j = surf%j(m)
    1896                 k = surf%k(m)
    1897 !
    1898 !--             If 2-m level is below the first grid level, MOST is
    1899 !--             used for calculation of 2-m temperature.
    1900                 IF ( surf%z_mo(m) > 2.0_wp )  THEN
    1901                    surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa     &
    1902                                    * ( LOG( 2.0_wp /  surf%z0h(m) )            &
    1903                                      - psi_h( 2.0_wp / surf%ol(m) )            &
    1904                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
    1905 !
    1906 !--             If 2-m level is above the first grid level, 2-m temperature
    1907 !--             is linearly interpolated between the two nearest vertical grid
    1908 !--             levels. Note, since 2-m temperature is only computed for
    1909 !--             horizontal upward-facing surfaces, only a vertical
    1910 !--             interpolation is necessary.
    1911                 ELSE
    1912 !
    1913 !--                zw(k-1) defines the height of the surface.
    1914                    kk = k
    1915                    DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp  .AND.  kk <= nzt )
    1916                       kk = kk + 1
    1917                    ENDDO               
    1918 !
    1919 !--                kk defines the index of the first grid level >= 2m.
    1920                    surf%pt_2m(m) = pt(kk-1,j,i) +                              &
    1921                                  ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *         &
    1922                                  ( pt(kk,j,i)       - pt(kk-1,j,i) ) /         &
    1923                                  ( zu(kk)           - zu(kk-1)     )
    1924                 ENDIF
    1925 
    1926              ENDDO
    1927          
    1928        
     1890
     1891             ENDDO
     1892
    19291893       END SELECT
    19301894
  • palm/trunk/SOURCE/surface_mod.f90

    r4329 r4331  
    2626! -----------------
    2727! $Id$
     28! -pt_2m - array is moved to diagnostic_output_quantities
     29!
     30! 4329 2019-12-10 15:46:36Z motisi
    2831! Renamed wall_flags_0 to wall_flags_static_0
    2932!
     
    289292       
    290293       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm             !< near surface air potential temperature at distance 10 cm from the surface (K)
    291        REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_2m               !< near surface air potential temperature at distance 2 m from the surface (K)     
    292294       
    293295       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  alpha_vg          !< coef. of Van Genuchten
     
    12811283!--    Salinity surface flux
    12821284       IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
    1283 !
    1284 !--    2-m potential temperature (for output quantity theta_2m*)
    1285        IF ( do_output_at_2m )  DEALLOCATE ( surfaces%pt_2m )
    12861285
    12871286    END SUBROUTINE deallocate_surface_attributes_h
     
    14121411!--    Salinity surface flux
    14131412       IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
    1414 !
    1415 !--    2-m potential temperature (for output quantity theta_2m*)
    1416        IF ( do_output_at_2m )  THEN
    1417           ALLOCATE ( surfaces%pt_2m(1:surfaces%ns) )
    1418           surfaces%pt_2m = -9999.0_wp  !< output array (for theta_2m*) must be initialized here,
    1419                                        !< otherwise simulation crash at do2d_at_begin with spinup=.F.
    1420        ENDIF
    14211413
    14221414    END SUBROUTINE allocate_surface_attributes_h
  • palm/trunk/SOURCE/time_integration_spinup.f90

    r4227 r4331  
    2525! -----------------
    2626! $Id$
     27! Enable output of diagnostic quantities, e.g. 2-m temperature
     28!
     29! 4227 2019-09-10 18:04:34Z gronemeier
    2730! implement new palm_date_time_mod
    2831!
     
    8184    USE cpulog,                                                                &
    8285        ONLY:  cpu_log, log_point_s
     86
     87    USE diagnostic_output_quantities_mod,                                      &
     88        ONLY:  doq_calculate
    8389
    8490    USE indices,                                                               &
     
    468474!--       2d-data output (cross-sections)
    469475          IF ( time_do2d_xy >= dt_do2d_xy )  THEN
     476             CALL doq_calculate
    470477             CALL data_output_2d( 'xy', 0 )
    471478             time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
     
    475482!--       3d-data output (volume data)
    476483          IF ( time_do3d >= dt_do3d )  THEN
     484             CALL doq_calculate
    477485             CALL data_output_3d( 0 )
    478486             time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
     
    560568
    561569       USE basic_constants_and_equations_mod,                                  &
    562        ONLY:  pi
     570           ONLY:  pi
    563571     
    564572       USE kinds
  • palm/trunk/SOURCE/write_restart_data_mod.f90

    r4301 r4331  
    2525! -----------------
    2626! $Id$
     27! - Move 2-m potential temperature to diagnostic_output_quantities.
     28!
     29! 4301 2019-11-22 12:09:09Z oliver.maas
    2730! removed recycling_yshift
    2831!
     
    942945       ENDIF
    943946
    944        IF ( ALLOCATED( pt_2m_av ) )  THEN
    945           CALL wrd_write_string( 'pt_2m_av' )
    946           WRITE ( 14 )  pt_2m_av
    947        ENDIF
    948 
    949947       IF ( humidity )  THEN
    950948
Note: See TracChangeset for help on using the changeset viewer.