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