phys/module_diag_nwp.F

   1 #if ( NMM_CORE == 1)
   2 MODULE module_diag_nwp
   3 CONTAINS
   4    SUBROUTINE diag_nwp_stub
   5    END SUBROUTINE diag_nwp_stub
   6 END MODULE module_diag_nwp
   7 #else
   8 !WRF:MEDIATION_LAYER:PHYSICS
   9 !
  10
  11 MODULE module_diag_nwp
  12       PRIVATE :: WGAMMA
  13       PRIVATE :: GAMMLN
  14 CONTAINS
  15    SUBROUTINE diagnostic_output_nwp(                                  &
  16                       config_flags,                                   &
  17                       ids,ide, jds,jde, kds,kde,                      &
  18                       ims,ime, jms,jme, kms,kme,                      &
  19                       ips,ipe, jps,jpe, kps,kpe,                      & ! patch  dims
  20                       i_start,i_end,j_start,j_end,kts,kte,num_tiles   &
  21                      ,u,v,temp,p8w                                    &
  22                      ,dt,xtime,sbw                                    &
  23                      ,mphysics_opt                                    &
  24                      ,gsfcgce_hail, gsfcgce_2ice                      &
  25                      ,mpuse_hail                                      &
  26                      ,nssl_cnoh, nssl_cnohl                           &
  27                      ,nssl_rho_qh, nssl_rho_qhl                       &
  28                      ,nssl_alphah, nssl_alphahl                       &
  29                      ,curr_secs2                                      &
  30                      ,nwp_diagnostics, diagflag                       &
  31                      ,history_interval                                &
  32                      ,itimestep                                       &
  33                      ,u10,v10,w                                       &
  34                      ,wspd10max                                       &
  35                      ,up_heli_max                                     &
  36                      ,w_up_max,w_dn_max                               &
  37                      ,znw,w_colmean                                   &
  38                      ,numcolpts,w_mean                                &
  39                      ,grpl_max,grpl_colint,refd_max,refl_10cm         &
  40                      ,hail_maxk1,hail_max2d                           &
  41                      ,qg_curr                                         &
  42                      ,ng_curr,qvolg_curr,qh_curr,nh_curr,qr_curr,nr_curr         & !  Optional (gthompsn)
  43                      ,rho,ph,phb,g                                    &
  44                      ,max_time_step,adaptive_ts                       &
  45                                                                       )
  46 !----------------------------------------------------------------------
  47
  48   USE module_configure, ONLY : grid_config_rec_type
  49
  50   USE module_state_description, ONLY :                                  &
  51       KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME, &
  52       WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, THOMPSONGH,         &
  53       MORR_TWO_MOMENT, GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME,           &
  54       MILBRANDT2MOM , CAMMGMPSCHEME, FULL_KHAIN_LYNN, MORR_TM_AERO,     &
  55       NSSL_2MOM, FAST_KHAIN_LYNN_SHPUND  !,MILBRANDT3MOM
  56   USE MODULE_MP_THOMPSON, ONLY: idx_bg1
  57
  58
  59    IMPLICIT NONE
  60 !======================================================================
  61 ! Definitions
  62 !-----------
  63 !-- DIAG_PRINT    print control: 0 - no diagnostics; 1 - dmudt only; 2 - all
  64 !-- DT            time step (second)
  65 !-- XTIME         forecast time
  66 !-- SBW           specified boundary width - used later
  67 !
  68 !-- P8W           3D pressure array at full eta levels
  69 !-- U             u component of wind - to be used later to compute k.e.
  70 !-- V             v component of wind - to be used later to compute k.e.
  71 !-- CURR_SECS2    Time (s) since the beginning of the restart
  72 !-- NWP_DIAGNOSTICS  = 1, compute hourly maximum fields
  73 !-- DIAGFLAG      logical flag to indicate if this is a history output time
  74 !-- U10, V10      10 m wind components
  75 !-- WSPD10MAX     10 m max wind speed
  76 !-- UP_HELI_MAX   max updraft helicity
  77 !-- W_UP_MAX      max updraft vertical velocity
  78 !-- W_DN_MAX      max downdraft vertical velocity
  79 !-- W_COLMEAN     column mean vertical velocity
  80 !-- NUMCOLPTS     no of column points
  81 !-- GRPL_MAX      max column-integrated graupel
  82 !-- GRPL_COLINT   column-integrated graupel
  83 !-- REF_MAX       max derived radar reflectivity
  84 !-- REFL_10CM     model computed 3D reflectivity
  85 !
  86 !-- ids           start index for i in domain
  87 !-- ide           end index for i in domain
  88 !-- jds           start index for j in domain
  89 !-- jde           end index for j in domain
  90 !-- kds           start index for k in domain
  91 !-- kde           end index for k in domain
  92 !-- ims           start index for i in memory
  93 !-- ime           end index for i in memory
  94 !-- jms           start index for j in memory
  95 !-- jme           end index for j in memory
  96 !-- ips           start index for i in patch
  97 !-- ipe           end index for i in patch
  98 !-- jps           start index for j in patch
  99 !-- jpe           end index for j in patch
 100 !-- kms           start index for k in memory
 101 !-- kme           end index for k in memory
 102 !-- i_start       start indices for i in tile
 103 !-- i_end         end indices for i in tile
 104 !-- j_start       start indices for j in tile
 105 !-- j_end         end indices for j in tile
 106 !-- kts           start index for k in tile
 107 !-- kte           end index for k in tile
 108 !-- num_tiles     number of tiles
 109 !
 110 !======================================================================
 111
 112       !  We are not changing any of the namelist settings.
 113
 114    TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags
 115
 116    INTEGER,      INTENT(IN   )    ::                             &
 117                                       ids,ide, jds,jde, kds,kde, &
 118                                       ims,ime, jms,jme, kms,kme, &
 119                                       ips,ipe, jps,jpe, kps,kpe, &
 120                                                         kts,kte, &
 121                                                       num_tiles
 122
 123    INTEGER, DIMENSION(num_tiles), INTENT(IN) ::                  &
 124      &           i_start,i_end,j_start,j_end
 125
 126    INTEGER,   INTENT(IN   )    ::   mphysics_opt
 127    INTEGER, INTENT(IN) :: gsfcgce_hail, gsfcgce_2ice, mpuse_hail
 128    REAL, INTENT(IN)    ::   nssl_cnoh, nssl_cnohl                &
 129                            ,nssl_rho_qh, nssl_rho_qhl            &
 130                            ,nssl_alphah, nssl_alphahl
 131    INTEGER,   INTENT(IN   )    ::   nwp_diagnostics
 132    LOGICAL,   INTENT(IN   )    ::   diagflag
 133
 134    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                 &
 135          INTENT(IN ) ::                                       u  &
 136                                                     ,         v  &
 137                                                     ,       p8w
 138
 139    REAL,  INTENT(IN   ) :: DT, XTIME
 140    INTEGER,  INTENT(IN   ) :: SBW
 141
 142    REAL, OPTIONAL, INTENT(IN)::  CURR_SECS2
 143
 144    INTEGER :: i,j,k,its,ite,jts,jte,ij
 145    INTEGER :: idp,jdp,irc,jrc,irnc,jrnc,isnh,jsnh
 146
 147    REAL              :: no_points
 148    REAL              :: dpsdt_sum, dmudt_sum, dardt_sum, drcdt_sum, drndt_sum
 149    REAL              :: hfx_sum, lh_sum, sfcevp_sum, rainc_sum, rainnc_sum, raint_sum
 150    REAL              :: dmumax, raincmax, rainncmax, snowhmax
 151    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
 152    CHARACTER*256     :: outstring
 153    CHARACTER*6       :: grid_str
 154
 155    INTEGER, INTENT(IN) ::                                        &
 156                                      history_interval,itimestep
 157
 158    REAL, DIMENSION( kms:kme ), INTENT(IN) ::                     &
 159                                                             znw
 160
 161    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) ::   &
 162                                                               w  &
 163                                                           ,temp  &
 164                                                        ,qg_curr  &
 165                                                            ,rho  &
 166                                                      ,refl_10cm  &
 167                                                         ,ph,phb
 168
 169    REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) ::    &
 170                            ng_curr, qvolg_curr, qh_curr, nh_curr        &
 171                                                ,qr_curr, nr_curr
 172
 173    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) ::            &
 174                                                             u10  &
 175                                                            ,v10
 176
 177    REAL, INTENT(IN) :: g
 178
 179    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) ::        &
 180                                                       wspd10max  &
 181                                                    ,up_heli_max  &
 182                                              ,w_up_max,w_dn_max  &
 183                                     ,w_colmean,numcolpts,w_mean  &
 184                                           ,grpl_max,grpl_colint  &
 185                                          ,hail_maxk1,hail_max2d  &
 186                                                       ,refd_max
 187
 188    REAL, DIMENSION(ims:ime,kms:kme,jms:jme)::                    &
 189                                       temp_qg, temp_ng, temp_qb  &
 190                                               ,temp_qr, temp_nr
 191    INTEGER, DIMENSION(ims:ime,kms:kme,jms:jme):: idx_bg
 192
 193    INTEGER :: idump
 194
 195    REAL :: wind_vel
 196    REAL :: depth
 197
 198       DOUBLE PRECISION:: hail_max
 199       REAL:: hail_max_sp
 200       DOUBLE PRECISION, PARAMETER:: thresh_conc = 0.0005d0                 ! number conc. of graupel/hail per cubic meter
 201       LOGICAL:: scheme_has_graupel
 202       INTEGER, PARAMETER:: ngbins=50
 203       DOUBLE PRECISION, DIMENSION(ngbins+1):: xxDx
 204       DOUBLE PRECISION, DIMENSION(ngbins):: xxDg, xdtg
 205       REAL:: xrho_g, xam_g, xxam_g, xbm_g, xmu_g
 206       REAL, DIMENSION(9), PARAMETER:: xrho_gh = (/ 50,100,200,300,400,500,600,700,800 /)
 207       REAL, DIMENSION(3):: cge, cgg
 208       DOUBLE PRECISION:: f_d, sum_ng, sum_t, lamg, ilamg, N0_g, lam_exp, N0exp
 209       DOUBLE PRECISION:: lamr, N0min
 210       REAL:: mvd_r, xslw1, ygra1, zans1
 211       INTEGER:: ng, n, k_0
 212
 213     REAL                                       :: time_from_output
 214     INTEGER                                    :: max_time_step
 215     LOGICAL                                    :: adaptive_ts
 216     LOGICAL                                    :: reset_arrays = .FALSE.
 217
 218 !-----------------------------------------------------------------
 219
 220     idump = (history_interval * 60.) / dt
 221
 222 ! IF ( MOD(itimestep, idump) .eq. 0 ) THEN
 223 !    WRITE(outstring,*) 'Computing PH0 for this domain with curr_secs2 = ', curr_secs2
 224 !    CALL wrf_message ( TRIM(outstring) )
 225
 226    time_from_output = mod( xtime, REAL(history_interval) )
 227
 228 !   print *,' history_interval     = ', history_interval
 229 !   print *,' itimestep            = ', itimestep
 230 !   print *,' idump                = ', idump
 231 !   print *,' xtime                = ', xtime
 232 !   print *,' time_from_output     = ', time_from_output
 233 !   print *,' max_time_step        = ', max_time_step
 234
 235    IF ( adaptive_ts .EQV. .TRUE. ) THEN
 236 !     if timestep is adaptive, use new resetting method;
 237 !     also, we are multiplying max_time_step with 1.05; because of rounding error,
 238 !     the time_from_output can become slightly larger than max_time_step when
 239 !     adaptive time step is maximized, in which case the if condition below fails to detect
 240 !     that we just wrote an output
 241        IF ( ( time_from_output .GT. 0 ) .AND. ( time_from_output .LE. ( ( max_time_step * 1.05 ) / 60. ) ) )  THEN
 242           reset_arrays = .TRUE.
 243        ENDIF
 244    ELSE
 245 !     if timestep is fixed, use original resetting method
 246       IF ( MOD((itimestep - 1), idump) .eq. 0 )  THEN
 247         reset_arrays = .TRUE.
 248       ENDIF
 249    ENDIF
 250
 251 !   print *,' reset_arrays = ', reset_arrays
 252    IF ( reset_arrays .EQV. .TRUE. ) THEN
 253
 254 !   print *,' Resetting NWP max arrays '
 255
 256 !  IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN
 257      WRITE(outstring,*) 'NSSL Diagnostics: Resetting max arrays for domain with dt = ', dt
 258      CALL wrf_debug ( 10,TRIM(outstring) )
 259
 260 !  !$OMP PARALLEL DO   &
 261 !  !$OMP PRIVATE ( ij )
 262      DO ij = 1 , num_tiles
 263        DO j=j_start(ij),j_end(ij)
 264        DO i=i_start(ij),i_end(ij)
 265          wspd10max(i,j)   = 0.
 266          up_heli_max(i,j) = 0.
 267          w_up_max(i,j)    = 0.
 268          w_dn_max(i,j)    = 0.
 269          w_mean(i,j)      = 0.
 270          grpl_max(i,j)    = 0.
 271          refd_max(i,j)    = 0.
 272          hail_maxk1(i,j)  = 0.
 273          hail_max2d(i,j)  = 0.
 274        ENDDO
 275        ENDDO
 276      ENDDO
 277 !  !$OMP END PARALLEL DO
 278      reset_arrays = .FALSE.
 279    ENDIF
 280
 281 !  !$OMP PARALLEL DO   &
 282 !  !$OMP PRIVATE ( ij )
 283    DO ij = 1 , num_tiles
 284      DO j=j_start(ij),j_end(ij)
 285      DO i=i_start(ij),i_end(ij)
 286
 287 ! Zero some accounting arrays that will be used below
 288
 289        w_colmean(i,j)   = 0.
 290        numcolpts(i,j)   = 0.
 291        grpl_colint(i,j) = 0.
 292      ENDDO
 293      ENDDO
 294    ENDDO
 295 !  !$OMP END PARALLEL DO
 296
 297 !  !$OMP PARALLEL DO   &
 298 !  !$OMP PRIVATE ( ij )
 299    DO ij = 1 , num_tiles
 300      DO j=j_start(ij),j_end(ij)
 301      DO k=kms,kme
 302      DO i=i_start(ij),i_end(ij)
 303
 304 ! Find vertical velocity max (up and down) below 400 mb
 305
 306        IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .GT. w_up_max(i,j) ) THEN
 307          w_up_max(i,j) = w(i,k,j)
 308        ENDIF
 309
 310        IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .LT. w_dn_max(i,j) ) THEN
 311          w_dn_max(i,j) = w(i,k,j)
 312        ENDIF
 313
 314 ! For the column mean vertical velocity calculation, first
 315 ! total the vertical velocity between sigma levels 0.5 and 0.8
 316
 317        IF ( znw(k) .GE. 0.5 .AND. znw(k) .LE. 0.8 ) THEN
 318          w_colmean(i,j) = w_colmean(i,j) + w(i,k,j)
 319          numcolpts(i,j) = numcolpts(i,j) + 1
 320        ENDIF
 321      ENDDO
 322      ENDDO
 323      ENDDO
 324    ENDDO
 325 !  !$OMP END PARALLEL DO
 326
 327 !  !$OMP PARALLEL DO   &
 328 !  !$OMP PRIVATE ( ij )
 329    DO ij = 1 , num_tiles
 330      DO j=j_start(ij),j_end(ij)
 331      DO k=kms,kme-1
 332      DO i=i_start(ij),i_end(ij)
 333
 334 ! Calculate the column integrated graupel
 335
 336        depth = ( ( ph(i,k+1,j) + phb(i,k+1,j) ) / g ) - &
 337                ( ( ph(i,k  ,j) + phb(i,k  ,j) ) / g )
 338        grpl_colint(i,j) = grpl_colint(i,j) + qg_curr(i,k,j) * depth * rho(i,k,j)
 339      ENDDO
 340      ENDDO
 341      ENDDO
 342    ENDDO
 343 !  !$OMP END PARALLEL DO
 344
 345 !  !$OMP PARALLEL DO   &
 346 !  !$OMP PRIVATE ( ij )
 347    DO ij = 1 , num_tiles
 348      DO j=j_start(ij),j_end(ij)
 349      DO i=i_start(ij),i_end(ij)
 350
 351 ! Calculate the max 10 m wind speed between output times
 352
 353        wind_vel = sqrt ( u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j) )
 354        IF ( wind_vel .GT. wspd10max(i,j) ) THEN
 355          wspd10max(i,j) = wind_vel
 356        ENDIF
 357
 358 ! Calculate the column mean vertical velocity between output times
 359
 360        w_mean(i,j) = w_mean(i,j) + w_colmean(i,j) / numcolpts(i,j)
 361
 362        IF ( MOD(itimestep, idump) .eq. 0 ) THEN
 363          w_mean(i,j) = w_mean(i,j) / idump
 364        ENDIF
 365
 366 ! Calculate the max column integrated graupel between output times
 367
 368        IF ( grpl_colint(i,j) .gt. grpl_max(i,j) ) THEN
 369           grpl_max(i,j) = grpl_colint(i,j)
 370        ENDIF
 371
 372    ! Calculate the max radar reflectivity between output times
 373
 374        IF ( refl_10cm(i,kms,j) .GT. refd_max(i,j) ) THEN
 375          refd_max(i,j) = refl_10cm(i,kms,j)
 376        ENDIF
 377      ENDDO
 378      ENDDO
 379    ENDDO
 380 !  !$OMP END PARALLEL DO
 381
 382 !+---+-----------------------------------------------------------------+
 383 !+---+-----------------------------------------------------------------+
 384 !..Calculate a maximum hail diameter from the characteristics of the
 385 !.. graupel category mixing ratio and number concentration (or hail, if
 386 !.. available).  This diagnostic uses the actual spectral distribution
 387 !.. assumptions, calculated by breaking the distribution into 50 bins
 388 !.. from 0.5mm to 7.5cm.  Once a minimum number concentration of thresh_conc (5e-4)
 389 !.. particle per cubic meter of air is reached, from the upper size
 390 !.. limit, then this bin is considered the max size.
 391 !+---+-----------------------------------------------------------------+
 392
 393    WRITE(outstring,*) 'GT-Diagnostics, computing max-hail diameter'
 394    CALL wrf_debug (100, TRIM(outstring))
 395
 396    IF (PRESENT(qh_curr)) THEN
 397    WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail mixing ratio'
 398    CALL wrf_debug (150, TRIM(outstring))
 399 !  !$OMP PARALLEL DO   &
 400 !  !$OMP PRIVATE ( ij )
 401    DO ij = 1 , num_tiles
 402      DO j=j_start(ij),j_end(ij)
 403      DO k=kms,kme-1
 404      DO i=i_start(ij),i_end(ij)
 405         temp_qg(i,k,j) = MAX(1.E-12, qh_curr(i,k,j)*rho(i,k,j))
 406      ENDDO
 407      ENDDO
 408      ENDDO
 409    ENDDO
 410 !  !$OMP END PARALLEL DO
 411    ELSE
 412 !  !$OMP PARALLEL DO   &
 413 !  !$OMP PRIVATE ( ij )
 414    DO ij = 1 , num_tiles
 415      DO j=j_start(ij),j_end(ij)
 416      DO k=kms,kme-1
 417      DO i=i_start(ij),i_end(ij)
 418         temp_qg(i,k,j) = MAX(1.E-12, qg_curr(i,k,j)*rho(i,k,j))
 419      ENDDO
 420      ENDDO
 421      ENDDO
 422    ENDDO
 423 !  !$OMP END PARALLEL DO
 424    ENDIF
 425
 426    IF (PRESENT(nh_curr)) THEN
 427    WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail number concentration'
 428    CALL wrf_debug (150, TRIM(outstring))
 429 !  !$OMP PARALLEL DO   &
 430 !  !$OMP PRIVATE ( ij )
 431    DO ij = 1 , num_tiles
 432      DO j=j_start(ij),j_end(ij)
 433      DO k=kms,kme-1
 434      DO i=i_start(ij),i_end(ij)
 435         temp_ng(i,k,j) = MAX(1.E-8, nh_curr(i,k,j)*rho(i,k,j))
 436      ENDDO
 437      ENDDO
 438      ENDDO
 439    ENDDO
 440 !  !$OMP END PARALLEL DO
 441    ELSEIF (PRESENT(ng_curr)) THEN
 442    WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel number concentration'
 443 !  !$OMP PARALLEL DO   &
 444 !  !$OMP PRIVATE ( ij )
 445    DO ij = 1 , num_tiles
 446      DO j=j_start(ij),j_end(ij)
 447      DO k=kms,kme-1
 448      DO i=i_start(ij),i_end(ij)
 449         temp_ng(i,k,j) = MAX(1.E-8, ng_curr(i,k,j)*rho(i,k,j))
 450      ENDDO
 451      ENDDO
 452      ENDDO
 453    ENDDO
 454 !  !$OMP END PARALLEL DO
 455    ELSE
 456 !  !$OMP PARALLEL DO   &
 457 !  !$OMP PRIVATE ( ij )
 458    DO ij = 1 , num_tiles
 459      DO j=j_start(ij),j_end(ij)
 460      DO k=kms,kme-1
 461      DO i=i_start(ij),i_end(ij)
 462         temp_ng(i,k,j) = 1.E-8
 463      ENDDO
 464      ENDDO
 465      ENDDO
 466    ENDDO
 467 !  !$OMP END PARALLEL DO
 468    ENDIF
 469
 470    IF (PRESENT(qvolg_curr)) THEN
 471    WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel volume mixing ratio'
 472    CALL wrf_debug (150, TRIM(outstring))
 473 !  !$OMP PARALLEL DO   &
 474 !  !$OMP PRIVATE ( ij )
 475    DO ij = 1 , num_tiles
 476      DO j=j_start(ij),j_end(ij)
 477      DO k=kms,kme-1
 478      DO i=i_start(ij),i_end(ij)
 479         temp_qb(i,k,j) = MAX(temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(9), qvolg_curr(i,k,j))
 480         temp_qb(i,k,j) = MIN(temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(1), temp_qb(i,k,j))
 481      ENDDO
 482      ENDDO
 483      ENDDO
 484    ENDDO
 485 !  !$OMP END PARALLEL DO
 486    ELSE
 487 !  !$OMP PARALLEL DO   &
 488 !  !$OMP PRIVATE ( ij )
 489    DO ij = 1 , num_tiles
 490      DO j=j_start(ij),j_end(ij)
 491      DO k=kms,kme-1
 492      DO i=i_start(ij),i_end(ij)
 493         idx_bg(i,k,j) = idx_bg1
 494         temp_qb(i,k,j) = temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(idx_bg(i,k,j))
 495      ENDDO
 496      ENDDO
 497      ENDDO
 498    ENDDO
 499 !  !$OMP END PARALLEL DO
 500    ENDIF
 501
 502
 503       ! Calculate the max hail size from graupel/hail parameters in microphysics scheme (gthompsn 12Mar2015)
 504       ! First, we do know multiple schemes have assumed inverse-exponential distribution with constant
 505       ! intercept parameter and particle density.  Please leave next 4 settings alone for common
 506       ! use and copy these and cge constants to re-assign per scheme if needed (e.g. NSSL).
 507
 508       xrho_g = 500.
 509       xam_g = 3.1415926536/6.0*xrho_g     ! Assumed m(D) = a*D**b, where a=PI/6*rho_g and b=3
 510       xbm_g = 3.                          ! in other words, spherical graupel/hail
 511       xmu_g = 0.
 512       scheme_has_graupel = .false.
 513
 514       !..Some constants. These *MUST* get changed below per scheme
 515       !.. *IF* either xbm_g or xmu_g value is changed from 3 and zero, respectively.
 516
 517       cge(1) = xbm_g + 1.
 518       cge(2) = xmu_g + 1.
 519       cge(3) = xbm_g + xmu_g + 1.
 520       do n = 1, 3
 521          cgg(n) = WGAMMA(cge(n))
 522       enddo
 523
 524    mp_select: SELECT CASE(mphysics_opt)
 525
 526      CASE (KESSLERSCHEME)
 527 !      nothing to do here
 528
 529      CASE (LINSCHEME)
 530        scheme_has_graupel = .true.
 531        xrho_g = 917.
 532        xam_g = 3.1415926536/6.0*xrho_g
 533        N0exp = 4.e4
 534 !      !$OMP PARALLEL DO   &
 535 !      !$OMP PRIVATE ( ij )
 536        DO ij = 1 , num_tiles
 537          DO j=j_start(ij),j_end(ij)
 538          DO k=kme-1, kms, -1
 539          DO i=i_start(ij),i_end(ij)
 540            if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
 541            lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
 542            temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2))
 543          ENDDO
 544          ENDDO
 545          ENDDO
 546        ENDDO
 547 !      !$OMP END PARALLEL DO
 548
 549      CASE (WSM3SCHEME)
 550 !      nothing to do here
 551
 552      CASE (WSM5SCHEME)
 553 !      nothing to do here
 554
 555      CASE (WSM6SCHEME)
 556        scheme_has_graupel = .true.
 557        xrho_g = 500.
 558        xam_g = 3.1415926536/6.0*xrho_g
 559        N0exp = 4.e6
 560 !      !$OMP PARALLEL DO   &
 561 !      !$OMP PRIVATE ( ij )
 562        DO ij = 1 , num_tiles
 563          DO j=j_start(ij),j_end(ij)
 564          DO k=kme-1, kms, -1
 565          DO i=i_start(ij),i_end(ij)
 566            if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
 567            lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
 568            temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2))
 569          ENDDO
 570          ENDDO
 571          ENDDO
 572        ENDDO
 573 !      !$OMP END PARALLEL DO
 574
 575      CASE (WDM5SCHEME)
 576 !      nothing to do here
 577
 578      CASE (WDM6SCHEME)
 579        scheme_has_graupel = .true.
 580        xrho_g = 500.
 581        N0exp = 4.e6
 582        if (mpuse_hail .eq. 1) then
 583          xrho_g = 700.
 584          N0exp = 4.e4
 585        endif
 586        xam_g = 3.1415926536/6.0*xrho_g
 587 !      !$OMP PARALLEL DO   &
 588 !      !$OMP PRIVATE ( ij )
 589        DO ij = 1 , num_tiles
 590          DO j=j_start(ij),j_end(ij)
 591          DO k=kme-1, kms, -1
 592          DO i=i_start(ij),i_end(ij)
 593            if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
 594            lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
 595            temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2))
 596          ENDDO
 597          ENDDO
 598          ENDDO
 599        ENDDO
 600 !      !$OMP END PARALLEL DO
 601
 602      CASE (GSFCGCESCHEME)
 603       if (gsfcgce_2ice.eq.0 .OR. gsfcgce_2ice.eq.2) then
 604        scheme_has_graupel = .true.
 605        if (gsfcgce_hail .eq. 1) then
 606           xrho_g = 900.
 607        else
 608           xrho_g = 400.
 609        endif
 610        xam_g = 3.1415926536/6.0*xrho_g
 611        N0exp = 4.e4
 612 !      !$OMP PARALLEL DO   &
 613 !      !$OMP PRIVATE ( ij )
 614        DO ij = 1 , num_tiles
 615          DO j=j_start(ij),j_end(ij)
 616          DO k=kme-1, kms, -1
 617          DO i=i_start(ij),i_end(ij)
 618            if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
 619            lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
 620            temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2))
 621          ENDDO
 622          ENDDO
 623          ENDDO
 624        ENDDO
 625 !      !$OMP END PARALLEL DO
 626       endif
 627
 628      CASE (SBU_YLINSCHEME)
 629 !      scheme_has_graupel = .true.  ! Can be calculated from rime fraction variable.
 630 !      no time to implement
 631
 632      CASE (ETAMPNEW)
 633 !      scheme_has_graupel = .true.  ! Can be calculated from rime fraction variable.
 634 !      no time to implement
 635
 636      CASE (THOMPSON, THOMPSONAERO)
 637
 638        scheme_has_graupel = .true.
 639        xmu_g = 0.
 640        cge(1) = xbm_g + 1.
 641        cge(2) = xmu_g + 1.
 642        cge(3) = xbm_g + xmu_g + 1.
 643        do n = 1, 3
 644           cgg(n) = WGAMMA(cge(n))
 645        enddo
 646
 647       if (.not.(PRESENT(ng_curr))) then
 648 !  !$OMP PARALLEL DO   &
 649 !  !$OMP PRIVATE ( ij )
 650       DO ij = 1 , num_tiles
 651        DO j=j_start(ij),j_end(ij)
 652        DO k=kms,kme-1
 653        DO i=i_start(ij),i_end(ij)
 654           temp_qr(i,k,j) = MAX(1.E-10, qr_curr(i,k,j)*rho(i,k,j))
 655           temp_nr(i,k,j) = MAX(1.E-8, nr_curr(i,k,j)*rho(i,k,j))
 656        ENDDO
 657        ENDDO
 658        ENDDO
 659       ENDDO
 660 !  !$OMP END PARALLEL DO
 661
 662 !  !$OMP PARALLEL DO   &
 663 !  !$OMP PRIVATE ( ij )
 664       DO ij = 1 , num_tiles
 665        DO j=j_start(ij),j_end(ij)
 666        DO i=i_start(ij),i_end(ij)
 667         DO k=kme-1, kms, -1
 668          ygra1 = alog10(max(1.E-9, temp_qg(i,k,j)))
 669          zans1 = 3.0 + 2./7.*(ygra1+8.)
 670          zans1 = MAX(2., MIN(zans1, 6.))
 671          N0exp = 10.**(zans1)
 672          lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
 673          lamg = lam_exp * (cgg(3)/cgg(2)/cgg(1))**(1./xbm_g)
 674          N0_g = N0exp/(cgg(2)*lam_exp) * lamg**cge(2)
 675          temp_ng(i,k,j) = MAX(1.E-8, N0_g*cgg(2)*lamg**(-cge(2)))
 676         ENDDO
 677        ENDDO
 678        ENDDO
 679       ENDDO
 680 !  !$OMP END PARALLEL DO
 681
 682       endif
 683
 684
 685      CASE (THOMPSONGH)
 686
 687        scheme_has_graupel = .true.
 688
 689 !  !$OMP PARALLEL DO   &
 690 !  !$OMP PRIVATE ( ij )
 691       DO ij = 1 , num_tiles
 692        DO j=j_start(ij),j_end(ij)
 693         DO k=kme-1, kms, -1
 694          DO i=i_start(ij),i_end(ij)
 695           idx_bg(i,k,j) = NINT(temp_qg(i,k,j)/temp_qb(i,k,j) *0.01)+1
 696           idx_bg(i,k,j) = MAX(1, MIN(idx_bg(i,k,j), 9))
 697         ENDDO
 698        ENDDO
 699        ENDDO
 700       ENDDO
 701 !  !$OMP END PARALLEL DO
 702
 703
 704 !    CASE (MORR_MILB_P3)
 705 !      scheme_has_graupel = .true.
 706 !      either Hugh or Jason need to implement.
 707
 708      CASE (MORR_TWO_MOMENT, MORR_TM_AERO)
 709        scheme_has_graupel = .true.
 710        xrho_g = 400.
 711        if (mpuse_hail .eq. 1) xrho_g = 900.
 712        xam_g = 3.1415926536/6.0*xrho_g
 713
 714      CASE (MILBRANDT2MOM)
 715        WRITE(outstring,*) 'GT-Debug, using Milbrandt2mom, which has 2-moment hail'
 716        CALL wrf_debug (150, TRIM(outstring))
 717        scheme_has_graupel = .true.
 718        xrho_g = 900.
 719        xam_g = 3.1415926536/6.0*xrho_g
 720
 721 !    CASE (MILBRANDT3MOM)
 722 !      coming in future?
 723
 724      CASE (NSSL_2MOM)
 725 ! Only treat 1-moment option here. 2- and 3-moment are now done in the microphysics
 726 !
 727        if ( config_flags%nssl_2moment_on == 0 ) then ! single-moment
 728
 729        scheme_has_graupel = .true.
 730        xrho_g = nssl_rho_qh
 731        N0exp = nssl_cnoh
 732        if (config_flags%nssl_hail_on==1) then
 733           xrho_g = nssl_rho_qhl
 734           N0exp = nssl_cnohl
 735        endif
 736        xam_g = 3.1415926536/6.0*xrho_g
 737
 738
 739        IF (config_flags%nssl_hail_on==1) THEN
 740          xmu_g = nssl_alphahl
 741        ELSE
 742          xmu_g = nssl_alphah
 743        ENDIF
 744
 745        if (xmu_g .NE. 0.) then
 746           cge(1) = xbm_g + 1.
 747           cge(2) = xmu_g + 1.
 748           cge(3) = xbm_g + xmu_g + 1.
 749           do n = 1, 3
 750              cgg(n) = WGAMMA(cge(n))
 751           enddo
 752        endif
 753
 754        ENDIF
 755
 756        ! NSSL scheme has many options, but, if single-moment, just fill
 757        ! in the number array for graupel from built-in assumptions.
 758
 759 !       if (.NOT.(PRESENT(nh_curr).OR.PRESENT(ng_curr)) ) then
 760        if ( config_flags%nssl_2moment_on == 0 ) then ! single-moment
 761 !      !$OMP PARALLEL DO   &
 762 !      !$OMP PRIVATE ( ij )
 763        DO ij = 1 , num_tiles
 764          DO j=j_start(ij),j_end(ij)
 765          DO k=kme-1, kms, -1
 766          DO i=i_start(ij),i_end(ij)
 767            if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
 768            lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
 769            temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2))
 770          ENDDO
 771          ENDDO
 772          ENDDO
 773        ENDDO
 774 !      !$OMP END PARALLEL DO
 775        endif
 776
 777 !    CASE (NSSL_3MOM)
 778 !      coming in future?
 779
 780      CASE (CAMMGMPSCHEME)
 781 !      nothing to do here
 782
 783      CASE (FULL_KHAIN_LYNN)
 784 !      explicit bin scheme so code below not applicable
 785 !      scheme authors need to implement if desired.
 786
 787      CASE (FAST_KHAIN_LYNN_SHPUND)
 788 !      explicit bin scheme so code below not applicable
 789 !      scheme authors need to implement if desired.
 790
 791      CASE DEFAULT
 792
 793    END SELECT mp_select
 794
 795
 796    if (scheme_has_graupel) then
 797
 798 !..Create bins of graupel/hail (from 500 microns up to 7.5 cm).
 799       xxDx(1) = 500.D-6
 800       xxDx(ngbins+1) = 0.075d0
 801       do n = 2, ngbins
 802          xxDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(ngbins) &
 803                   *DLOG(xxDx(ngbins+1)/xxDx(1)) +DLOG(xxDx(1)))
 804       enddo
 805       do n = 1, ngbins
 806          xxDg(n) = DSQRT(xxDx(n)*xxDx(n+1))
 807          xdtg(n) = xxDx(n+1) - xxDx(n)
 808       enddo
 809
 810
 811 !  !$OMP PARALLEL DO   &
 812 !  !$OMP PRIVATE ( ij )
 813       DO ij = 1 , num_tiles
 814         DO j=j_start(ij),j_end(ij)
 815         DO k=kms,kme-1
 816         DO i=i_start(ij),i_end(ij)
 817            if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
 818            IF (PRESENT(qvolg_curr)) THEN
 819             xxam_g = 3.1415926536/6.0*xrho_gh(idx_bg(i,k,j))
 820             IF (xrho_gh(idx_bg(i,k,j)).lt.350.0) CYCLE
 821            ELSE
 822             xxam_g = xam_g
 823            ENDIF
 824            lamg = (xxam_g*cgg(3)/cgg(2)*temp_ng(i,k,j)/temp_qg(i,k,j))**(1./xbm_g)
 825            N0_g = temp_ng(i,k,j)/cgg(2)*lamg**cge(2)
 826            sum_ng = 0.0d0
 827            sum_t  = 0.0d0
 828            do ng = ngbins, 1, -1
 829               f_d = N0_g*xxDg(ng)**xmu_g * DEXP(-lamg*xxDg(ng))*xdtg(ng)
 830               sum_ng = sum_ng + f_d
 831               if (sum_ng .gt. thresh_conc) then
 832                  exit
 833               endif
 834               sum_t = sum_ng
 835            enddo
 836            if (ng .ge. ngbins) then
 837               hail_max = xxDg(ngbins)
 838            elseif (xxDg(ng+1) .gt. 1.E-3) then
 839               hail_max = xxDg(ng) - (sum_ng-thresh_conc)/(sum_ng-sum_t) &
 840      &                              * (xxDg(ng)-xxDg(ng+1))
 841            else
 842               hail_max = 1.E-4
 843            endif
 844 !          if (hail_max .gt. 1E-2) then
 845 !           WRITE(outstring,*) 'GT-Debug-Hail, ', hail_max*1000.
 846 !           CALL wrf_debug (350, TRIM(outstring))
 847 !          endif
 848            hail_max_sp = hail_max
 849            if (k.eq.kms) then
 850               hail_maxk1(i,j) = MAX(hail_maxk1(i,j), hail_max_sp)
 851            endif
 852            hail_max2d(i,j) = MAX(hail_max2d(i,j), hail_max_sp)
 853         ENDDO
 854         ENDDO
 855         ENDDO
 856       ENDDO
 857 !  !$OMP END PARALLEL DO
 858
 859    endif
 860
 861    END SUBROUTINE diagnostic_output_nwp
 862
 863 !+---+-----------------------------------------------------------------+
 864       REAL FUNCTION GAMMLN(XX)
 865 !     --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
 866       IMPLICIT NONE
 867       REAL, INTENT(IN):: XX
 868       DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0
 869       DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
 870                COF = (/76.18009172947146D0, -86.50532032941677D0, &
 871                        24.01409824083091D0, -1.231739572450155D0, &
 872                       .1208650973866179D-2, -.5395239384953D-5/)
 873       DOUBLE PRECISION:: SER,TMP,X,Y
 874       INTEGER:: J
 875
 876       X=XX
 877       Y=X
 878       TMP=X+5.5D0
 879       TMP=(X+0.5D0)*LOG(TMP)-TMP
 880       SER=1.000000000190015D0
 881       DO 11 J=1,6
 882         Y=Y+1.D0
 883         SER=SER+COF(J)/Y
 884 11    CONTINUE
 885       GAMMLN=TMP+LOG(STP*SER/X)
 886       END FUNCTION GAMMLN
 887 !  (C) Copr. 1986-92 Numerical Recipes Software 2.02
 888 !+---+-----------------------------------------------------------------+
 889       REAL FUNCTION WGAMMA(y)
 890
 891       IMPLICIT NONE
 892       REAL, INTENT(IN):: y
 893
 894       WGAMMA = EXP(GAMMLN(y))
 895
 896       END FUNCTION WGAMMA
 897 !+---+-----------------------------------------------------------------+
 898
 899 END MODULE module_diag_nwp
 900 #endif