#if ( NMM_CORE == 1) MODULE module_diag_nwp CONTAINS SUBROUTINE diag_nwp_stub END SUBROUTINE diag_nwp_stub END MODULE module_diag_nwp #else !WRF:MEDIATION_LAYER:PHYSICS ! MODULE module_diag_nwp PRIVATE :: WGAMMA PRIVATE :: GAMMLN CONTAINS SUBROUTINE diagnostic_output_nwp( & config_flags, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & ! patch dims i_start,i_end,j_start,j_end,kts,kte,num_tiles & ,u,v,temp,p8w & ,dt,xtime,sbw & ,mphysics_opt & ,gsfcgce_hail, gsfcgce_2ice & ,mpuse_hail & ,nssl_cnoh, nssl_cnohl & ,nssl_rho_qh, nssl_rho_qhl & ,nssl_alphah, nssl_alphahl & ,curr_secs2 & ,nwp_diagnostics, diagflag & ,history_interval & ,itimestep & ,u10,v10,w & ,wspd10max & ,up_heli_max & ,w_up_max,w_dn_max & ,znw,w_colmean & ,numcolpts,w_mean & ,grpl_max,grpl_colint,refd_max,refl_10cm & ,hail_maxk1,hail_max2d & ,qg_curr & ,ng_curr,qvolg_curr,qh_curr,nh_curr,qr_curr,nr_curr & ! Optional (gthompsn) ,rho,ph,phb,g & ,max_time_step,adaptive_ts & ) !---------------------------------------------------------------------- USE module_configure, ONLY : grid_config_rec_type USE module_state_description, ONLY : & KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME, & WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, THOMPSONGH, & MORR_TWO_MOMENT, GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, & MILBRANDT2MOM , CAMMGMPSCHEME, FULL_KHAIN_LYNN, MORR_TM_AERO, & NSSL_2MOM, FAST_KHAIN_LYNN_SHPUND !,MILBRANDT3MOM USE MODULE_MP_THOMPSON, ONLY: idx_bg1 IMPLICIT NONE !====================================================================== ! Definitions !----------- !-- DIAG_PRINT print control: 0 - no diagnostics; 1 - dmudt only; 2 - all !-- DT time step (second) !-- XTIME forecast time !-- SBW specified boundary width - used later ! !-- P8W 3D pressure array at full eta levels !-- U u component of wind - to be used later to compute k.e. !-- V v component of wind - to be used later to compute k.e. !-- CURR_SECS2 Time (s) since the beginning of the restart !-- NWP_DIAGNOSTICS = 1, compute hourly maximum fields !-- DIAGFLAG logical flag to indicate if this is a history output time !-- U10, V10 10 m wind components !-- WSPD10MAX 10 m max wind speed !-- UP_HELI_MAX max updraft helicity !-- W_UP_MAX max updraft vertical velocity !-- W_DN_MAX max downdraft vertical velocity !-- W_COLMEAN column mean vertical velocity !-- NUMCOLPTS no of column points !-- GRPL_MAX max column-integrated graupel !-- GRPL_COLINT column-integrated graupel !-- REF_MAX max derived radar reflectivity !-- REFL_10CM model computed 3D reflectivity ! !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- ips start index for i in patch !-- ipe end index for i in patch !-- jps start index for j in patch !-- jpe end index for j in patch !-- kms start index for k in memory !-- kme end index for k in memory !-- i_start start indices for i in tile !-- i_end end indices for i in tile !-- j_start start indices for j in tile !-- j_end end indices for j in tile !-- kts start index for k in tile !-- kte end index for k in tile !-- num_tiles number of tiles ! !====================================================================== ! We are not changing any of the namelist settings. TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags INTEGER, INTENT(IN ) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & kts,kte, & num_tiles INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & & i_start,i_end,j_start,j_end INTEGER, INTENT(IN ) :: mphysics_opt INTEGER, INTENT(IN) :: gsfcgce_hail, gsfcgce_2ice, mpuse_hail REAL, INTENT(IN) :: nssl_cnoh, nssl_cnohl & ,nssl_rho_qh, nssl_rho_qhl & ,nssl_alphah, nssl_alphahl INTEGER, INTENT(IN ) :: nwp_diagnostics LOGICAL, INTENT(IN ) :: diagflag REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(IN ) :: u & , v & , p8w REAL, INTENT(IN ) :: DT, XTIME INTEGER, INTENT(IN ) :: SBW REAL, OPTIONAL, INTENT(IN):: CURR_SECS2 INTEGER :: i,j,k,its,ite,jts,jte,ij INTEGER :: idp,jdp,irc,jrc,irnc,jrnc,isnh,jsnh REAL :: no_points REAL :: dpsdt_sum, dmudt_sum, dardt_sum, drcdt_sum, drndt_sum REAL :: hfx_sum, lh_sum, sfcevp_sum, rainc_sum, rainnc_sum, raint_sum REAL :: dmumax, raincmax, rainncmax, snowhmax LOGICAL, EXTERNAL :: wrf_dm_on_monitor CHARACTER*256 :: outstring CHARACTER*6 :: grid_str INTEGER, INTENT(IN) :: & history_interval,itimestep REAL, DIMENSION( kms:kme ), INTENT(IN) :: & znw REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: & w & ,temp & ,qg_curr & ,rho & ,refl_10cm & ,ph,phb REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: & ng_curr, qvolg_curr, qh_curr, nh_curr & ,qr_curr, nr_curr REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: & u10 & ,v10 REAL, INTENT(IN) :: g REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: & wspd10max & ,up_heli_max & ,w_up_max,w_dn_max & ,w_colmean,numcolpts,w_mean & ,grpl_max,grpl_colint & ,hail_maxk1,hail_max2d & ,refd_max REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: & temp_qg, temp_ng, temp_qb & ,temp_qr, temp_nr INTEGER, DIMENSION(ims:ime,kms:kme,jms:jme):: idx_bg INTEGER :: idump REAL :: wind_vel REAL :: depth DOUBLE PRECISION:: hail_max REAL:: hail_max_sp DOUBLE PRECISION, PARAMETER:: thresh_conc = 0.0005d0 ! number conc. of graupel/hail per cubic meter LOGICAL:: scheme_has_graupel INTEGER, PARAMETER:: ngbins=50 DOUBLE PRECISION, DIMENSION(ngbins+1):: xxDx DOUBLE PRECISION, DIMENSION(ngbins):: xxDg, xdtg REAL:: xrho_g, xam_g, xxam_g, xbm_g, xmu_g REAL, DIMENSION(9), PARAMETER:: xrho_gh = (/ 50,100,200,300,400,500,600,700,800 /) REAL, DIMENSION(3):: cge, cgg DOUBLE PRECISION:: f_d, sum_ng, sum_t, lamg, ilamg, N0_g, lam_exp, N0exp DOUBLE PRECISION:: lamr, N0min REAL:: mvd_r, xslw1, ygra1, zans1 INTEGER:: ng, n, k_0 REAL :: time_from_output INTEGER :: max_time_step LOGICAL :: adaptive_ts LOGICAL :: reset_arrays = .FALSE. !----------------------------------------------------------------- idump = (history_interval * 60.) / dt ! IF ( MOD(itimestep, idump) .eq. 0 ) THEN ! WRITE(outstring,*) 'Computing PH0 for this domain with curr_secs2 = ', curr_secs2 ! CALL wrf_message ( TRIM(outstring) ) time_from_output = mod( xtime, REAL(history_interval) ) ! print *,' history_interval = ', history_interval ! print *,' itimestep = ', itimestep ! print *,' idump = ', idump ! print *,' xtime = ', xtime ! print *,' time_from_output = ', time_from_output ! print *,' max_time_step = ', max_time_step IF ( adaptive_ts .EQV. .TRUE. ) THEN ! if timestep is adaptive, use new resetting method; ! also, we are multiplying max_time_step with 1.05; because of rounding error, ! the time_from_output can become slightly larger than max_time_step when ! adaptive time step is maximized, in which case the if condition below fails to detect ! that we just wrote an output IF ( ( time_from_output .GT. 0 ) .AND. ( time_from_output .LE. ( ( max_time_step * 1.05 ) / 60. ) ) ) THEN reset_arrays = .TRUE. ENDIF ELSE ! if timestep is fixed, use original resetting method IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN reset_arrays = .TRUE. ENDIF ENDIF ! print *,' reset_arrays = ', reset_arrays IF ( reset_arrays .EQV. .TRUE. ) THEN ! print *,' Resetting NWP max arrays ' ! IF ( MOD((itimestep - 1), idump) .eq. 0 ) THEN WRITE(outstring,*) 'NSSL Diagnostics: Resetting max arrays for domain with dt = ', dt CALL wrf_debug ( 10,TRIM(outstring) ) ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) wspd10max(i,j) = 0. up_heli_max(i,j) = 0. w_up_max(i,j) = 0. w_dn_max(i,j) = 0. w_mean(i,j) = 0. grpl_max(i,j) = 0. refd_max(i,j) = 0. hail_maxk1(i,j) = 0. hail_max2d(i,j) = 0. ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO reset_arrays = .FALSE. ENDIF ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) ! Zero some accounting arrays that will be used below w_colmean(i,j) = 0. numcolpts(i,j) = 0. grpl_colint(i,j) = 0. ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme DO i=i_start(ij),i_end(ij) ! Find vertical velocity max (up and down) below 400 mb IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .GT. w_up_max(i,j) ) THEN w_up_max(i,j) = w(i,k,j) ENDIF IF ( p8w(i,k,j) .GT. 40000. .AND. w(i,k,j) .LT. w_dn_max(i,j) ) THEN w_dn_max(i,j) = w(i,k,j) ENDIF ! For the column mean vertical velocity calculation, first ! total the vertical velocity between sigma levels 0.5 and 0.8 IF ( znw(k) .GE. 0.5 .AND. znw(k) .LE. 0.8 ) THEN w_colmean(i,j) = w_colmean(i,j) + w(i,k,j) numcolpts(i,j) = numcolpts(i,j) + 1 ENDIF ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) ! Calculate the column integrated graupel depth = ( ( ph(i,k+1,j) + phb(i,k+1,j) ) / g ) - & ( ( ph(i,k ,j) + phb(i,k ,j) ) / g ) grpl_colint(i,j) = grpl_colint(i,j) + qg_curr(i,k,j) * depth * rho(i,k,j) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) ! Calculate the max 10 m wind speed between output times wind_vel = sqrt ( u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j) ) IF ( wind_vel .GT. wspd10max(i,j) ) THEN wspd10max(i,j) = wind_vel ENDIF ! Calculate the column mean vertical velocity between output times w_mean(i,j) = w_mean(i,j) + w_colmean(i,j) / numcolpts(i,j) IF ( MOD(itimestep, idump) .eq. 0 ) THEN w_mean(i,j) = w_mean(i,j) / idump ENDIF ! Calculate the max column integrated graupel between output times IF ( grpl_colint(i,j) .gt. grpl_max(i,j) ) THEN grpl_max(i,j) = grpl_colint(i,j) ENDIF ! Calculate the max radar reflectivity between output times IF ( refl_10cm(i,kms,j) .GT. refd_max(i,j) ) THEN refd_max(i,j) = refl_10cm(i,kms,j) ENDIF ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ !..Calculate a maximum hail diameter from the characteristics of the !.. graupel category mixing ratio and number concentration (or hail, if !.. available). This diagnostic uses the actual spectral distribution !.. assumptions, calculated by breaking the distribution into 50 bins !.. from 0.5mm to 7.5cm. Once a minimum number concentration of thresh_conc (5e-4) !.. particle per cubic meter of air is reached, from the upper size !.. limit, then this bin is considered the max size. !+---+-----------------------------------------------------------------+ WRITE(outstring,*) 'GT-Diagnostics, computing max-hail diameter' CALL wrf_debug (100, TRIM(outstring)) IF (PRESENT(qh_curr)) THEN WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail mixing ratio' CALL wrf_debug (150, TRIM(outstring)) ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) temp_qg(i,k,j) = MAX(1.E-12, qh_curr(i,k,j)*rho(i,k,j)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ELSE ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) temp_qg(i,k,j) = MAX(1.E-12, qg_curr(i,k,j)*rho(i,k,j)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ENDIF IF (PRESENT(nh_curr)) THEN WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has hail number concentration' CALL wrf_debug (150, TRIM(outstring)) ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) temp_ng(i,k,j) = MAX(1.E-8, nh_curr(i,k,j)*rho(i,k,j)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ELSEIF (PRESENT(ng_curr)) THEN WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel number concentration' ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) temp_ng(i,k,j) = MAX(1.E-8, ng_curr(i,k,j)*rho(i,k,j)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ELSE ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) temp_ng(i,k,j) = 1.E-8 ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ENDIF IF (PRESENT(qvolg_curr)) THEN WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel volume mixing ratio' CALL wrf_debug (150, TRIM(outstring)) ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) temp_qb(i,k,j) = MAX(temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(9), qvolg_curr(i,k,j)) temp_qb(i,k,j) = MIN(temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(1), temp_qb(i,k,j)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ELSE ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) idx_bg(i,k,j) = idx_bg1 temp_qb(i,k,j) = temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(idx_bg(i,k,j)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ENDIF ! Calculate the max hail size from graupel/hail parameters in microphysics scheme (gthompsn 12Mar2015) ! First, we do know multiple schemes have assumed inverse-exponential distribution with constant ! intercept parameter and particle density. Please leave next 4 settings alone for common ! use and copy these and cge constants to re-assign per scheme if needed (e.g. NSSL). xrho_g = 500. xam_g = 3.1415926536/6.0*xrho_g ! Assumed m(D) = a*D**b, where a=PI/6*rho_g and b=3 xbm_g = 3. ! in other words, spherical graupel/hail xmu_g = 0. scheme_has_graupel = .false. !..Some constants. These *MUST* get changed below per scheme !.. *IF* either xbm_g or xmu_g value is changed from 3 and zero, respectively. cge(1) = xbm_g + 1. cge(2) = xmu_g + 1. cge(3) = xbm_g + xmu_g + 1. do n = 1, 3 cgg(n) = WGAMMA(cge(n)) enddo mp_select: SELECT CASE(mphysics_opt) CASE (KESSLERSCHEME) ! nothing to do here CASE (LINSCHEME) scheme_has_graupel = .true. xrho_g = 917. xam_g = 3.1415926536/6.0*xrho_g N0exp = 4.e4 ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kme-1, kms, -1 DO i=i_start(ij),i_end(ij) if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO CASE (WSM3SCHEME) ! nothing to do here CASE (WSM5SCHEME) ! nothing to do here CASE (WSM6SCHEME) scheme_has_graupel = .true. xrho_g = 500. xam_g = 3.1415926536/6.0*xrho_g N0exp = 4.e6 ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kme-1, kms, -1 DO i=i_start(ij),i_end(ij) if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO CASE (WDM5SCHEME) ! nothing to do here CASE (WDM6SCHEME) scheme_has_graupel = .true. xrho_g = 500. N0exp = 4.e6 if (mpuse_hail .eq. 1) then xrho_g = 700. N0exp = 4.e4 endif xam_g = 3.1415926536/6.0*xrho_g ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kme-1, kms, -1 DO i=i_start(ij),i_end(ij) if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO CASE (GSFCGCESCHEME) if (gsfcgce_2ice.eq.0 .OR. gsfcgce_2ice.eq.2) then scheme_has_graupel = .true. if (gsfcgce_hail .eq. 1) then xrho_g = 900. else xrho_g = 400. endif xam_g = 3.1415926536/6.0*xrho_g N0exp = 4.e4 ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kme-1, kms, -1 DO i=i_start(ij),i_end(ij) if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO endif CASE (SBU_YLINSCHEME) ! scheme_has_graupel = .true. ! Can be calculated from rime fraction variable. ! no time to implement CASE (ETAMPNEW) ! scheme_has_graupel = .true. ! Can be calculated from rime fraction variable. ! no time to implement CASE (THOMPSON, THOMPSONAERO) scheme_has_graupel = .true. xmu_g = 0. cge(1) = xbm_g + 1. cge(2) = xmu_g + 1. cge(3) = xbm_g + xmu_g + 1. do n = 1, 3 cgg(n) = WGAMMA(cge(n)) enddo if (.not.(PRESENT(ng_curr))) then ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) temp_qr(i,k,j) = MAX(1.E-10, qr_curr(i,k,j)*rho(i,k,j)) temp_nr(i,k,j) = MAX(1.E-8, nr_curr(i,k,j)*rho(i,k,j)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO i=i_start(ij),i_end(ij) DO k=kme-1, kms, -1 ygra1 = alog10(max(1.E-9, temp_qg(i,k,j))) zans1 = 3.0 + 2./7.*(ygra1+8.) zans1 = MAX(2., MIN(zans1, 6.)) N0exp = 10.**(zans1) lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) lamg = lam_exp * (cgg(3)/cgg(2)/cgg(1))**(1./xbm_g) N0_g = N0exp/(cgg(2)*lam_exp) * lamg**cge(2) temp_ng(i,k,j) = MAX(1.E-8, N0_g*cgg(2)*lamg**(-cge(2))) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO endif CASE (THOMPSONGH) scheme_has_graupel = .true. ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kme-1, kms, -1 DO i=i_start(ij),i_end(ij) idx_bg(i,k,j) = NINT(temp_qg(i,k,j)/temp_qb(i,k,j) *0.01)+1 idx_bg(i,k,j) = MAX(1, MIN(idx_bg(i,k,j), 9)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO ! CASE (MORR_MILB_P3) ! scheme_has_graupel = .true. ! either Hugh or Jason need to implement. CASE (MORR_TWO_MOMENT, MORR_TM_AERO) scheme_has_graupel = .true. xrho_g = 400. if (mpuse_hail .eq. 1) xrho_g = 900. xam_g = 3.1415926536/6.0*xrho_g CASE (MILBRANDT2MOM) WRITE(outstring,*) 'GT-Debug, using Milbrandt2mom, which has 2-moment hail' CALL wrf_debug (150, TRIM(outstring)) scheme_has_graupel = .true. xrho_g = 900. xam_g = 3.1415926536/6.0*xrho_g ! CASE (MILBRANDT3MOM) ! coming in future? CASE (NSSL_2MOM) ! Only treat 1-moment option here. 2- and 3-moment are now done in the microphysics ! if ( config_flags%nssl_2moment_on == 0 ) then ! single-moment scheme_has_graupel = .true. xrho_g = nssl_rho_qh N0exp = nssl_cnoh if (config_flags%nssl_hail_on==1) then xrho_g = nssl_rho_qhl N0exp = nssl_cnohl endif xam_g = 3.1415926536/6.0*xrho_g IF (config_flags%nssl_hail_on==1) THEN xmu_g = nssl_alphahl ELSE xmu_g = nssl_alphah ENDIF if (xmu_g .NE. 0.) then cge(1) = xbm_g + 1. cge(2) = xmu_g + 1. cge(3) = xbm_g + xmu_g + 1. do n = 1, 3 cgg(n) = WGAMMA(cge(n)) enddo endif ENDIF ! NSSL scheme has many options, but, if single-moment, just fill ! in the number array for graupel from built-in assumptions. ! if (.NOT.(PRESENT(nh_curr).OR.PRESENT(ng_curr)) ) then if ( config_flags%nssl_2moment_on == 0 ) then ! single-moment ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kme-1, kms, -1 DO i=i_start(ij),i_end(ij) if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1)) temp_ng(i,k,j) = N0exp*cgg(2)*lam_exp**(-cge(2)) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO endif ! CASE (NSSL_3MOM) ! coming in future? CASE (CAMMGMPSCHEME) ! nothing to do here CASE (FULL_KHAIN_LYNN) ! explicit bin scheme so code below not applicable ! scheme authors need to implement if desired. CASE (FAST_KHAIN_LYNN_SHPUND) ! explicit bin scheme so code below not applicable ! scheme authors need to implement if desired. CASE DEFAULT END SELECT mp_select if (scheme_has_graupel) then !..Create bins of graupel/hail (from 500 microns up to 7.5 cm). xxDx(1) = 500.D-6 xxDx(ngbins+1) = 0.075d0 do n = 2, ngbins xxDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(ngbins) & *DLOG(xxDx(ngbins+1)/xxDx(1)) +DLOG(xxDx(1))) enddo do n = 1, ngbins xxDg(n) = DSQRT(xxDx(n)*xxDx(n+1)) xdtg(n) = xxDx(n+1) - xxDx(n) enddo ! !$OMP PARALLEL DO & ! !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles DO j=j_start(ij),j_end(ij) DO k=kms,kme-1 DO i=i_start(ij),i_end(ij) if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE IF (PRESENT(qvolg_curr)) THEN xxam_g = 3.1415926536/6.0*xrho_gh(idx_bg(i,k,j)) IF (xrho_gh(idx_bg(i,k,j)).lt.350.0) CYCLE ELSE xxam_g = xam_g ENDIF lamg = (xxam_g*cgg(3)/cgg(2)*temp_ng(i,k,j)/temp_qg(i,k,j))**(1./xbm_g) N0_g = temp_ng(i,k,j)/cgg(2)*lamg**cge(2) sum_ng = 0.0d0 sum_t = 0.0d0 do ng = ngbins, 1, -1 f_d = N0_g*xxDg(ng)**xmu_g * DEXP(-lamg*xxDg(ng))*xdtg(ng) sum_ng = sum_ng + f_d if (sum_ng .gt. thresh_conc) then exit endif sum_t = sum_ng enddo if (ng .ge. ngbins) then hail_max = xxDg(ngbins) elseif (xxDg(ng+1) .gt. 1.E-3) then hail_max = xxDg(ng) - (sum_ng-thresh_conc)/(sum_ng-sum_t) & & * (xxDg(ng)-xxDg(ng+1)) else hail_max = 1.E-4 endif ! if (hail_max .gt. 1E-2) then ! WRITE(outstring,*) 'GT-Debug-Hail, ', hail_max*1000. ! CALL wrf_debug (350, TRIM(outstring)) ! endif hail_max_sp = hail_max if (k.eq.kms) then hail_maxk1(i,j) = MAX(hail_maxk1(i,j), hail_max_sp) endif hail_max2d(i,j) = MAX(hail_max2d(i,j), hail_max_sp) ENDDO ENDDO ENDDO ENDDO ! !$OMP END PARALLEL DO endif END SUBROUTINE diagnostic_output_nwp !+---+-----------------------------------------------------------------+ REAL FUNCTION GAMMLN(XX) ! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0. IMPLICIT NONE REAL, INTENT(IN):: XX DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: & COF = (/76.18009172947146D0, -86.50532032941677D0, & 24.01409824083091D0, -1.231739572450155D0, & .1208650973866179D-2, -.5395239384953D-5/) DOUBLE PRECISION:: SER,TMP,X,Y INTEGER:: J X=XX Y=X TMP=X+5.5D0 TMP=(X+0.5D0)*LOG(TMP)-TMP SER=1.000000000190015D0 DO 11 J=1,6 Y=Y+1.D0 SER=SER+COF(J)/Y 11 CONTINUE GAMMLN=TMP+LOG(STP*SER/X) END FUNCTION GAMMLN ! (C) Copr. 1986-92 Numerical Recipes Software 2.02 !+---+-----------------------------------------------------------------+ REAL FUNCTION WGAMMA(y) IMPLICIT NONE REAL, INTENT(IN):: y WGAMMA = EXP(GAMMLN(y)) END FUNCTION WGAMMA !+---+-----------------------------------------------------------------+ END MODULE module_diag_nwp #endif