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