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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.