diff --git a/src/physics/mp_driver.f90 b/src/physics/mp_driver.f90 index 462b81c8..06450430 100644 --- a/src/physics/mp_driver.f90 +++ b/src/physics/mp_driver.f90 @@ -3,12 +3,14 @@ module module_mp_driver use domain_interface, only: domain_t use module_mp_thompson, only: thompson_init, mp_gt_driver + private + public :: microphysics + logical :: initialized = .false. contains - subroutine mp_init(domain) + subroutine mp_init() implicit none - type(domain_t), intent(inout) :: domain call thompson_init() @@ -39,9 +41,8 @@ subroutine process_subdomain(domain, dt, its,ite, jts,jte, kts,kte) RAINNC=domain%accumulated_precipitation,& SNOWNC=domain%accumulated_snowfall, & has_reqc=0, has_reqi=0, has_reqs=0, & - ids=domain%ids,ide=domain%ide, & ! domain dims - jds=domain%jds,jde=domain%jde, & - kds=domain%kds,kde=domain%kde, & + ide=domain%ide, & ! domain dims + jde=domain%jde, & ims=domain%ims,ime=domain%ime, & ! memory dims jms=domain%jms,jme=domain%jme, & kms=domain%kms,kme=domain%kme, & @@ -98,7 +99,7 @@ subroutine microphysics(domain, dt, halo, subset) real, intent(in) :: dt integer, intent(in), optional :: halo, subset - if (.not.initialized) call mp_init(domain) + if (.not.initialized) call mp_init() if (present(subset)) then call process_subdomain(domain, dt, & @@ -130,9 +131,8 @@ subroutine microphysics(domain, dt, halo, subset) RAINNC=domain%accumulated_precipitation,& SNOWNC=domain%accumulated_snowfall, & has_reqc=0, has_reqi=0, has_reqs=0, & - ids=domain%ids,ide=domain%ide, & ! domain dims - jds=domain%jds,jde=domain%jde, & - kds=domain%kds,kde=domain%kde, & + ide=domain%ide, & ! domain dims + jde=domain%jde, & ims=domain%ims,ime=domain%ime, & ! memory dims jms=domain%jms,jme=domain%jme, & kms=domain%kms,kme=domain%kme, & diff --git a/src/physics/mp_thompson.f90 b/src/physics/mp_thompson.f90 index f3539a01..3cc357ea 100644 --- a/src/physics/mp_thompson.f90 +++ b/src/physics/mp_thompson.f90 @@ -47,19 +47,9 @@ MODULE module_mp_thompson use co_util, only : co_bcast use timer_interface, only : timer_t -! USE module_wrf_error -! USE module_mp_radar -! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) -! USE module_dm, ONLY : wrf_dm_max_real -! #endif - IMPLICIT NONE - - LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false. - LOGICAL, PRIVATE:: is_aerosol_aware = .false. - !$omp threadprivate(is_aerosol_aware) - LOGICAL, PARAMETER, PRIVATE:: dustyIce = .true. - LOGICAL, PARAMETER, PRIVATE:: homogIce = .true. + private + public :: mp_gt_driver, thompson_init INTEGER, PARAMETER, PRIVATE:: IFDRY = 0 REAL, PARAMETER, PRIVATE:: T_0 = 273.15 @@ -318,10 +308,6 @@ MODULE module_mp_thompson ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) REAL, DIMENSION(ntb_art), PARAMETER, PRIVATE:: & ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) - REAL, DIMENSION(ntb_arr), PARAMETER, PRIVATE:: & - ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) - REAL, DIMENSION(ntb_ark), PARAMETER, PRIVATE:: & - ta_Ka = (/0.2, 0.4, 0.6, 0.8/) !..Lookup tables for IN concentration (/m**3) from 0.001 to 1000/Liter. REAL, DIMENSION(ntb_IN), PARAMETER, PRIVATE:: & @@ -391,171 +377,17 @@ MODULE module_mp_thompson REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me -!+---+ -!+---+-----------------------------------------------------------------+ -!..END DECLARATIONS -!+---+-----------------------------------------------------------------+ -!+---+ -!ctrlL - CONTAINS - SUBROUTINE thompson_init(hgt, nwfa2d, nwfa, nifa, dx, dy, & - is_start, & - ids, ide, jds, jde, kds, kde, & - ims, ime, jms, jme, kms, kme, & - its, ite, jts, jte, kts, kte) + SUBROUTINE thompson_init() IMPLICIT NONE - INTEGER, OPTIONAL, INTENT(IN):: ids,ide, jds,jde, kds,kde, & - ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte - REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN):: hgt - -!..OPTIONAL variables that control application of aerosol-aware scheme - - REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: nwfa, nifa - REAL, DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: nwfa2d - REAL, OPTIONAL, INTENT(IN) :: DX, DY - LOGICAL, OPTIONAL, INTENT(IN) :: is_start - CHARACTER*256:: mp_debug - - INTEGER:: i, j, k, l, m, n - REAL:: h_01, niIN3, niCCN3, max_test - LOGICAL:: micro_init, has_CCN, has_IN + LOGICAL:: micro_init type(timer_t) :: timer - is_aerosol_aware = .FALSE. micro_init = .FALSE. - has_CCN = .FALSE. - has_IN = .FALSE. - - ! write(mp_debug,*) ' DEBUG checking column of hgt ', its+1,jts+1 - ! CALL wrf_debug(250, mp_debug) - ! do k = kts, kte - ! write(mp_debug,*) ' DEBUGT k, hgt = ', k, hgt(its+1,k,jts+1) - ! CALL wrf_debug(250, mp_debug) - ! enddo - - if (PRESENT(nwfa2d) .AND. PRESENT(nwfa) .AND. PRESENT(nifa)) is_aerosol_aware = .TRUE. - - if (is_aerosol_aware) then - -!..Check for existing aerosol data, both CCN and IN aerosols. If missing -!.. fill in just a basic vertical profile, somewhat boundary-layer following. - -! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) -! max_test = wrf_dm_max_real ( MAXVAL(nwfa(its:ite-1,:,jts:jte-1)) ) -! #else - max_test = MAXVAL ( nwfa(its:ite-1,:,jts:jte-1) ) -! #endif - - if (max_test .lt. eps) then - ! write(mp_debug,*) ' Apparently there are no initial CCN aerosols.' - ! CALL wrf_debug(100, mp_debug) - ! write(mp_debug,*) ' checked column at point (i,j) = ', its,jts - ! CALL wrf_debug(100, mp_debug) - do j = jts, min(jde-1,jte) - do i = its, min(ide-1,ite) - if (hgt(i,1,j).le.1000.0) then - h_01 = 0.8 - elseif (hgt(i,1,j).ge.2500.0) then - h_01 = 0.01 - else - h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) - endif - niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01 - nwfa(i,1,j) = naCCN1+naCCN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niCCN3) - do k = 2, kte - nwfa(i,k,j) = naCCN1+naCCN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niCCN3) - enddo - enddo - enddo - else - has_CCN = .TRUE. - ! write(mp_debug,*) ' Apparently initial CCN aerosols are present.' - ! CALL wrf_debug(100, mp_debug) - ! write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nwfa(its,:,jts)) - ! CALL wrf_debug(100, mp_debug) - endif - - -! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) -! max_test = wrf_dm_max_real ( MAXVAL(nifa(its:ite-1,:,jts:jte-1)) ) -! #else - max_test = MAXVAL ( nifa(its:ite-1,:,jts:jte-1) ) -! #endif - - if (max_test .lt. eps) then - ! write(mp_debug,*) ' Apparently there are no initial IN aerosols.' - ! CALL wrf_debug(100, mp_debug) - ! write(mp_debug,*) ' checked column at point (i,j) = ', its,jts - ! CALL wrf_debug(100, mp_debug) - do j = jts, min(jde-1,jte) - do i = its, min(ide-1,ite) - if (hgt(i,1,j).le.1000.0) then - h_01 = 0.8 - elseif (hgt(i,1,j).ge.2500.0) then - h_01 = 0.01 - else - h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) - endif - niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01 - nifa(i,1,j) = naIN1+naIN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niIN3) - do k = 2, kte - nifa(i,k,j) = naIN1+naIN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niIN3) - enddo - enddo - enddo - else - has_IN = .TRUE. - ! write(mp_debug,*) ' Apparently initial IN aerosols are present.' - ! CALL wrf_debug(100, mp_debug) - ! write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nifa(its,:,jts)) - ! CALL wrf_debug(100, mp_debug) - endif - -!..Capture initial state lowest level CCN aerosol data in 2D array. - -! do j = jts, min(jde-1,jte) -! do i = its, min(ide-1,ite) -! nwfa2d(i,j) = nwfa(i,kts,j) -! enddo -! enddo - -!..Scale the lowest level aerosol data into an emissions rate. This is -!.. very far from ideal, but need higher emissions where larger amount -!.. of existing and lesser emissions where not already lots of aerosols -!.. for first-order simplistic approach. Later, proper connection to -!.. emission inventory would be better, but, for now, scale like this: -!.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per kg per second -!.. Nwfa=500 per cc, emit 0.875E5 aerosols per kg per second -!.. Nwfa=5000 per cc, emit 0.875E6 aerosols per kg per second -!.. for a grid with 20km spacing and scale accordingly for other spacings. - - if (is_start) then - if (SQRT(DX*DY)/20000.0 .ge. 1.0) then - h_01 = 0.875 - else - h_01 = (0.875 + 0.125*((20000.-SQRT(DX*DY))/16000.)) * SQRT(DX*DY)/20000. - endif - ! write(mp_debug,*) ' aerosol surface flux emission scale factor is: ', h_01 - ! CALL wrf_debug(100, mp_debug) - do j = jts, min(jde-1,jte) - do i = its, min(ide-1,ite) - nwfa2d(i,j) = 10.0**(LOG10(nwfa(i,kts,j)*1.E-6)-3.69897) - nwfa2d(i,j) = nwfa2d(i,j)*h_01 * 1.E6 - enddo - enddo -! else -! write(mp_debug,*) ' sample (lower-left) aerosol surface flux emission rate: ', nwfa2d(1,1) -! CALL wrf_debug(100, mp_debug) - endif - - endif - !..Allocate space for lookup tables (J. Michalakes 2009Jun08). @@ -950,44 +782,16 @@ SUBROUTINE thompson_init(hgt, nwfa2d, nwfa, nifa, dx, dy, & enddo enddo - ! CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ') - ! WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') & - ! ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g - ! CALL wrf_debug(150, wrf_err_message) - -!..Read a static file containing CCN activation of aerosols. The -!.. data were created from a parcel model by Feingold & Heymsfield with -!.. further changes by Eidhammer and Kriedenweis. - if (is_aerosol_aware) then - ! CALL wrf_debug(200, ' calling table_ccnAct routine') - call table_ccnAct - endif - !..Collision efficiency between rain/snow and cloud water. - ! CALL wrf_debug(200, ' creating qc collision eff tables') call table_Efrw call table_Efsw !..Drop evaporation. - ! CALL wrf_debug(200, ' creating rain evap table') call table_dropEvap !..Initialize various constants for computing radar reflectivity. - ! xam_r = am_r - ! xbm_r = bm_r - ! xmu_r = mu_r - ! xam_s = am_s - ! xbm_s = bm_s - ! xmu_s = mu_s - ! xam_g = am_g - ! xbm_g = bm_g - ! xmu_g = mu_g - ! call radar_init - - if (.not. iiwarm) then !..Rain collecting graupel & graupel collecting rain. - ! CALL wrf_debug(200, ' creating rain collecting graupel table') call timer%reset() call timer%start() call qr_acr_qg @@ -998,7 +802,6 @@ SUBROUTINE thompson_init(hgt, nwfa2d, nwfa, nifa, dx, dy, & !..Rain collecting snow & snow collecting rain. - ! CALL wrf_debug(200, ' creating rain collecting snow table') call timer%reset() call timer%start() call qr_acr_qs @@ -1008,7 +811,6 @@ SUBROUTINE thompson_init(hgt, nwfa2d, nwfa, nifa, dx, dy, & endif !..Cloud water and rain freezing (Bigg, 1953). - ! CALL wrf_debug(200, ' creating freezing of water drops table') call timer%reset() call timer%start() call freezeH2O @@ -1018,7 +820,6 @@ SUBROUTINE thompson_init(hgt, nwfa2d, nwfa, nifa, dx, dy, & endif !..Conversion of some ice mass into snow category. - ! CALL wrf_debug(200, ' creating ice converting to snow table') call timer%reset() call timer%start() call qi_aut_qs @@ -1029,71 +830,41 @@ SUBROUTINE thompson_init(hgt, nwfa2d, nwfa, nifa, dx, dy, & endif - ! CALL wrf_debug(150, ' ... DONE microphysical lookup tables') - - endif - END SUBROUTINE thompson_init !+---+-----------------------------------------------------------------+ -!ctrlL -!+---+-----------------------------------------------------------------+ !..This is a wrapper routine designed to transfer values from 3D to 1D. !+---+-----------------------------------------------------------------+ - SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & - nwfa, nifa, nwfa2d, & - th, pii, p, w, dz, dt_in, itimestep, & - RAINNC, RAINNCV, & - SNOWNC, SNOWNCV, & - GRAUPELNC, GRAUPELNCV, SR, & -! #if ( WRF_CHEM == 1 ) -! rainprod, evapprod, & -! #endif - refl_10cm, diagflag, do_radar_ref, & - re_cloud, re_ice, re_snow, & + SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, & + th, pii, p, w, dz, dt_in, & + RAINNC, & + SNOWNC, & has_reqc, has_reqi, has_reqs, & - ids,ide, jds,jde, kds,kde, & ! domain dims + ide,jde, & ! domain dims ims,ime, jms,jme, kms,kme, & ! memory dims its,ite, jts,jte, kts,kte) ! tile dims implicit none !..Subroutine arguments - INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & + INTEGER, INTENT(IN):: ide,jde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & qv, qc, qr, qi, qs, qg, ni, nr, th - REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & - nc, nwfa, nifa - REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d - REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & - re_cloud, re_ice, re_snow INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs -! #if ( WRF_CHEM == 1 ) -! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & -! rainprod, evapprod -! #endif REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: & pii, p, w, dz REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & RAINNC REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: & - RAINNCV, SR, SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV - REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & - refl_10cm + SNOWNC REAL, INTENT(IN):: dt_in - INTEGER, optional, INTENT(IN):: itimestep -!..Local variables REAL, DIMENSION(kts:kte):: & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nr1d, nc1d, nwfa1d, nifa1d, & - t1d, p1d, w1d, dz1d, rho, dBZ + t1d, p1d, w1d, dz1d, rho REAL, DIMENSION(kts:kte):: re_qc1d, re_qi1d, re_qs1d -! #if ( WRF_CHEM == 1 ) -! REAL, DIMENSION(kts:kte):: & -! rainprod1d, evapprod1d -! #endif REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic REAL:: dt, pptrain, pptsnow, pptgraul, pptice REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max @@ -1103,10 +874,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr INTEGER:: i_start, j_start, i_end, j_end - LOGICAL, OPTIONAL, INTENT(IN) :: diagflag - INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref CHARACTER*256:: mp_debug - integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM !+---+ !$OMP PARALLEL DEFAULT(PRIVATE) FIRSTPRIVATE(ids,ide,jds,jde,kds,kde,ims,ime,jms,jme, & @@ -1121,15 +889,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & i_end = MIN(ite, ide-1) j_end = MIN(jte, jde-1) -!..For idealized testing by developer. -! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. & -! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then -! i_start = its + 2 -! i_end = ite -! j_start = jts -! j_end = jte -! endif - dt = dt_in qc_max = 0. @@ -1164,25 +923,14 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & mp_debug(i:i) = char(0) enddo - if (.NOT. is_aerosol_aware .AND. PRESENT(nc) .AND. PRESENT(nwfa) & - .AND. PRESENT(nifa) .AND. PRESENT(nwfa2d)) then - ! write(mp_debug,*) 'WARNING, nc-nwfa-nifa-nwfa2d present but is_aerosol_aware is FALSE' - ! CALL wrf_debug(0, mp_debug) - endif - !$omp do schedule(guided) j_loop: do j = j_start, j_end - ! print*, this_image(), j-j_start, omp_get_thread_num(), omp_get_num_threads() i_loop: do i = i_start, i_end pptrain = 0. pptsnow = 0. pptgraul = 0. pptice = 0. - IF ( PRESENT (RAINNCV) ) RAINNCV(i,j) = 0. - IF ( PRESENT (snowncv) ) SNOWNCV(i,j) = 0. - IF ( PRESENT (graupelncv) ) GRAUPELNCV(i,j) = 0. - IF ( PRESENT (SR) ) SR(i,j) = 0. do k = kts, kte t1d(k) = th(i,k,j)*pii(i,k,j) @@ -1198,63 +946,31 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ni1d(k) = ni(i,k,j) nr1d(k) = nr(i,k,j) enddo - if (is_aerosol_aware) then - do k = kts, kte - nc1d(k) = nc(i,k,j) - nwfa1d(k) = nwfa(i,k,j) - nifa1d(k) = nifa(i,k,j) - enddo - nwfa1 = nwfa2d(i,j) - else - do k = kts, kte - rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) - nc1d(k) = Nt_c/rho(k) - nwfa1d(k) = 11.1E6/rho(k) - nifa1d(k) = naIN1*0.01/rho(k) - enddo - nwfa1 = 11.1E6 - endif + do k = kts, kte + rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) + nc1d(k) = Nt_c/rho(k) + nwfa1d(k) = 11.1E6/rho(k) + nifa1d(k) = naIN1*0.01/rho(k) + enddo + nwfa1 = 11.1E6 call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, & pptrain, pptsnow, pptgraul, pptice, & -! #if ( WRF_CHEM == 1 ) -! rainprod1d, evapprod1d, & -! #endif - kts, kte, dt, i, j) + kts, kte, dt) pcp_ra(i,j) = pptrain pcp_sn(i,j) = pptsnow pcp_gr(i,j) = pptgraul pcp_ic(i,j) = pptice - IF ( PRESENT (RAINNCV) )RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice - IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN - SNOWNCV(i,j) = pptsnow + pptice + IF (PRESENT(snownc) ) THEN SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice ENDIF - IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN - GRAUPELNCV(i,j) = pptgraul - GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul - ENDIF - IF ( PRESENT(SR) )SR(i,j) = (pptsnow + pptgraul + pptice)/(pptrain + pptsnow + pptgraul + pptice+1.e-12) - - !..Reset lowest model level to initial state aerosols (fake sfc source). !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol !.. number tendency (number per kg per second). - if (is_aerosol_aware) then -!-GT nwfa1d(kts) = nwfa1 - nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt_in - - do k = kts, kte - nc(i,k,j) = nc1d(k) - nwfa(i,k,j) = nwfa1d(k) - nifa(i,k,j) = nifa1d(k) - enddo - endif - do k = kts, kte qv(i,k,j) = qv1d(k) qc(i,k,j) = qc1d(k) @@ -1265,87 +981,50 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ni(i,k,j) = ni1d(k) nr(i,k,j) = nr1d(k) th(i,k,j) = t1d(k)/pii(i,k,j) -! #if ( WRF_CHEM == 1 ) -! rainprod(i,k,j) = rainprod1d(k) -! evapprod(i,k,j) = evapprod1d(k) -! #endif if (qc1d(k) .gt. qc_max) then imax_qc = i jmax_qc = j kmax_qc = k qc_max = qc1d(k) - elseif (qc1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) endif if (qr1d(k) .gt. qr_max) then imax_qr = i jmax_qr = j kmax_qr = k qr_max = qr1d(k) - ! elseif (qr1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) endif if (nr1d(k) .gt. nr_max) then imax_nr = i jmax_nr = j kmax_nr = k nr_max = nr1d(k) - ! elseif (nr1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) endif if (qs1d(k) .gt. qs_max) then imax_qs = i jmax_qs = j kmax_qs = k qs_max = qs1d(k) - ! elseif (qs1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) endif if (qi1d(k) .gt. qi_max) then imax_qi = i jmax_qi = j kmax_qi = k qi_max = qi1d(k) - ! elseif (qi1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) endif if (qg1d(k) .gt. qg_max) then imax_qg = i jmax_qg = j kmax_qg = k qg_max = qg1d(k) - ! elseif (qg1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) endif if (ni1d(k) .gt. ni_max) then imax_ni = i jmax_ni = j kmax_ni = k ni_max = ni1d(k) - ! elseif (ni1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) endif if (qv1d(k) .lt. 0.0) then - ! write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), & - ! ' at i,j,k=', i,j,k - ! CALL wrf_debug(150, mp_debug) if (k.lt.kte-2 .and. k.gt.kts+1) then - ! write(mp_debug,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) - ! CALL wrf_debug(150, mp_debug) qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j))) else qv(i,k,j) = 1.E-7 @@ -1354,16 +1033,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & endif enddo - ! IF ( PRESENT (diagflag) ) THEN - ! if (diagflag .and. do_radar_ref == 1) then - ! call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & - ! t1d, p1d, dBZ, kts, kte, i, j) - ! do k = kts, kte - ! refl_10cm(i,k,j) = MAX(-35., dBZ(k)) - ! enddo - ! endif - ! ENDIF - IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN do k = kts, kte re_qc1d(k) = 2.49E-6 @@ -1372,11 +1041,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & enddo call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & re_qc1d, re_qi1d, re_qs1d, kts, kte) - do k = kts, kte - if (present(re_cloud)) re_cloud(i,k,j) = MAX(2.49E-6, MIN(re_qc1d(k), 50.E-6)) - if (present(re_ice)) re_ice(i,k,j) = MAX(4.99E-6, MIN(re_qi1d(k), 125.E-6)) - if (present(re_snow)) re_snow(i,k,j) = MAX(9.99E-6, MIN(re_qs1d(k), 999.E-6)) - enddo ENDIF enddo i_loop @@ -1384,27 +1048,12 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & !$omp end do !$omp end parallel -! DEBUG - GT - ! write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', & - ! 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', & - ! 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', & - ! 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', & - ! 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', & - ! 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', & - ! 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', & - ! 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')' - ! CALL wrf_debug(150, mp_debug) -! END DEBUG - GT - do i = 1, 256 mp_debug(i:i) = char(0) enddo END SUBROUTINE mp_gt_driver -!+---+-----------------------------------------------------------------+ -!ctrlL -!+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ !.. This subroutine computes the moisture tendencies of water vapor, !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel. @@ -1416,25 +1065,18 @@ END SUBROUTINE mp_gt_driver subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, & pptrain, pptsnow, pptgraul, pptice, & -! #if ( WRF_CHEM == 1 ) -! rainprod, evapprod, & -! #endif - kts, kte, dt, ii, jj) + kts, kte, dt) implicit none !..Sub arguments - INTEGER, INTENT(IN):: kts, kte, ii, jj + INTEGER, INTENT(IN):: kts, kte REAL, DIMENSION(kts:kte), INTENT(INOUT):: & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nr1d, nc1d, nwfa1d, nifa1d, t1d REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice REAL, INTENT(IN):: dt -! #if ( WRF_CHEM == 1 ) -! REAL, DIMENSION(kts:kte), INTENT(INOUT):: & -! rainprod, evapprod -! #endif !..Local variables REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, & @@ -1469,8 +1111,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & prg_gcw, prg_rci, prg_rcs, & prg_rcg, prg_ihm - DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0 - REAL, DIMENSION(kts:kte):: temp, pres, qv REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa REAL, DIMENSION(kts:kte):: rho, rhof, rhof2 @@ -1486,21 +1126,20 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c - REAL:: rgvm, delta_tp, orho, lfus2 + REAL:: delta_tp, orho, lfus2 REAL, DIMENSION(5):: onstep DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg DOUBLE PRECISION:: lami, ilami, ilamc - REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m - DOUBLE PRECISION:: Dr_star, Dc_star + REAL:: xDc, Dc_b, Dc_g, xDi, xDs, xDg REAL:: zeta1, zeta, taud, tau - REAL:: stoke_r, stoke_s, stoke_g, stoke_i + REAL:: stoke_g REAL:: vti, vtr, vts, vtg, vtc REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, & vtck, vtnck REAL, DIMENSION(kts:kte):: vts_boost REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow - REAL:: a_, b_, loga_, A1, A2, tf - REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat + REAL:: a_, b_, loga_, tf + REAL:: tempc, tc0 REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr REAL:: xsat, rate_max, sump, ratio REAL:: clap, fcd, dfcd @@ -1510,26 +1149,18 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & REAL:: Ef_ra, Ef_sa, Ef_ga REAL:: dtsave, odts, odt, odzq, hgt_agl REAL:: xslw1, ygra1, zans1, eva_factor - INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq + INTEGER:: k, n, nn, nstep, k_0, IT, iexfrq INTEGER, DIMENSION(5):: ksed1 INTEGER:: nir, nis, nig, nii, nic, niin INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, & - idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in + idx_i1, idx_i, idx_c, idx, idx_n, idx_in - LOGICAL:: melti, no_micro + LOGICAL:: no_micro LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg LOGICAL:: debug_flag - CHARACTER*256:: mp_debug INTEGER:: nu_c -!+---+ - debug_flag = .false. -! if (ii.eq.901 .and. jj.eq.379) debug_flag = .true. - ! if(debug_flag) then - ! write(mp_debug, *) 'DEBUG INFO, mp_thompson at (i,j) ', ii, ', ', jj - ! CALL wrf_debug(550, mp_debug) - ! endif no_micro = .true. dtsave = dt @@ -1634,12 +1265,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pnd_scd(k) = 0. pnd_gcd(k) = 0. enddo -! #if ( WRF_CHEM == 1 ) -! do k = kts, kte -! rainprod(k) = 0. -! evapprod(k) = 0. -! enddo -! #endif !..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments. do k = kts, kte @@ -1679,7 +1304,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) & / am_r*lamc**bm_r) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + nc(k) = Nt_c else qc1d(k) = 0.0 nc1d(k) = 0.0 @@ -1763,18 +1388,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif enddo -!+---+-----------------------------------------------------------------+ -! if (debug_flag) then -! write(mp_debug,*) 'DEBUG-VERBOSE at (i,j) ', ii, ', ', jj -! CALL wrf_debug(550, mp_debug) -! do k = kts, kte -! write(mp_debug, '(a,i3,f8.2,1x,f7.2,1x, 11(1x,e13.6))') & -! & 'VERBOSE: ', k, pres(k)*0.01, temp(k)-273.15, qv(k), rc(k), rr(k), ri(k), rs(k), rg(k), nc(k), nr(k), ni(k), nwfa(k), nifa(k) -! CALL wrf_debug(550, mp_debug) -! enddo -! endif -!+---+-----------------------------------------------------------------+ - !+---+-----------------------------------------------------------------+ !..Derive various thermodynamic variables frequently used. !.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from @@ -1821,7 +1434,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope, and useful moments for snow. !+---+-----------------------------------------------------------------+ - if (.not. iiwarm) then do k = kts, kte if (.not. L_qs(k)) CYCLE tc0 = MIN(-0.1, temp(k)-273.15) @@ -1932,8 +1544,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) enddo - endif - !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope values for rain. !+---+-----------------------------------------------------------------+ @@ -1953,10 +1563,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !..Rain self-collection follows Seifert, 1994 and drop break-up !.. follows Verlinde and Cotton, 1993. RAIN2M if (L_qr(k) .and. mvd_r(k).gt. D0r) then -!-GT Ef_rr = 1.0 -!-GT if (mvd_r(k) .gt. 1500.0E-6) then - Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6)) -!-GT endif + Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6)) pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k) endif @@ -2019,7 +1626,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !..Compute all frozen hydrometeor species' process terms. !+---+-----------------------------------------------------------------+ - if (.not. iiwarm) then do k = kts, kte vts_boost(k) = 1.5 @@ -2313,11 +1919,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !.. Implemented by T. Eidhammer and G. Thompson 2012Dec18 !+---+-----------------------------------------------------------------+ - if (dustyIce .AND. is_aerosol_aware) then - xni = iceDeMott(tempc,qvs(k),qvs(k),qvsi(k),rho(k),nifa(k)) - else - xni = 1.0 *1000. ! Default is 1.0 per Liter - endif + xni = 1.0 *1000. ! Default is 1.0 per Liter !..Ice nuclei lookup table index. if (xni.gt. Nt_IN(1)) then @@ -2362,11 +1964,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !.. we may need to relax the temperature and ssati constraints. if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps & .and. temp(k).lt.253.15) ) then - if (dustyIce .AND. is_aerosol_aware) then - xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k)) - else - xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k)))) - endif + xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k)))) xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts pri_inu(k) = MIN(DBLE(rate_max), xm0i*pni_inu(k)) @@ -2375,13 +1973,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !..Freezing of aqueous aerosols based on Koop et al (2001, Nature) xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave - if (is_aerosol_aware .AND. homogIce .AND. (xni.le.500.E3) & - & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then - xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave) - pni_iha(k) = xnc*odts - pri_iha(k) = MIN(DBLE(rate_max), xm0i*0.1*pni_iha(k)) - pni_iha(k) = pri_iha(k)/(xm0i*0.1) - endif !+---+------------------ END NEW ICE NUCLEATION -----------------------+ @@ -2529,12 +2120,9 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delQvs(k)) & * N0_g(k)*(t1_qg_me*ilamg(k)**cge(10) & + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11)) -!-GT prr_gml(k) = prr_gml(k) + 4218.*olfus*tempc & -!-GT * (prr_rcg(k)+prg_gcw(k)) prr_gml(k) = MIN(DBLE(rg(k)*odts), MAX(0.D0, prr_gml(k))) pnr_gml(k) = N0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) & ! RAIN2M * prr_gml(k) * 10.0**(-0.5*tempc) -! if (tempc.gt.7.5 .or. rg(k).lt.0.005E-3) pnr_gml(k)=0.0 if (ssati(k).lt. 0.) then prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs & @@ -2557,7 +2145,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif enddo - endif !+---+-----------------------------------------------------------------+ !..Ensure we do not deplete more hydrometeor species than exists. @@ -2668,19 +2255,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & orho = 1./rho(k) lfus2 = lsub - lvap(k) -!..Aerosol number tendency - if (is_aerosol_aware) then - nwfaten(k) = nwfaten(k) - (pna_rca(k) + pna_sca(k) & - + pna_gca(k) + pni_iha(k)) * orho - nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) & - + pnd_gcd(k)) * orho - if (dustyIce) then - nifaten(k) = nifaten(k) - pni_inu(k)*orho - else - nifaten(k) = 0. - endif - endif - !..Water vapor tendency qvten(k) = qvten(k) + (-pri_inu(k) - pri_iha(k) - pri_ide(k) & - prs_ide(k) - prs_sde(k) - prg_gde(k)) & @@ -2859,7 +2433,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & if ((qc1d(k) + qcten(k)*DT) .gt. R1) then rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k) nc(k) = MAX(2., (nc1d(k) + ncten(k)*DT)*rho(k)) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + nc(k) = Nt_c L_qc(k) = .true. else rc(k) = R1 @@ -2919,7 +2493,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !..With tendency-updated mixing ratios, recalculate snow moments and !.. intercepts/slopes of graupel and rain. !+---+-----------------------------------------------------------------+ - if (.not. iiwarm) then do k = kts, kte if (.not. L_qs(k)) CYCLE tc0 = MIN(-0.1, temp(k)-273.15) @@ -2997,8 +2570,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) enddo - endif - !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope values for rain. !+---+-----------------------------------------------------------------+ @@ -3033,62 +2604,9 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & prw_vcd(k) = clap*odt !+---+-----------------------------------------------------------------+ ! DROPLET NUCLEATION if (clap .gt. eps) then - if (is_aerosol_aware) then - xnc = MAX(2., activ_ncloud(temp(k), w1d(k), nwfa(k))) - else - xnc = Nt_c - endif + xnc = Nt_c pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho -!+---+-----------------------------------------------------------------+ ! EVAPORATION - elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.E-6 .AND. is_aerosol_aware) then - tempc = temp(k) - 273.15 - otemp = 1./temp(k) - rvs = rho(k)*qvs(k) - rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.) - rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) & - *otemp*(lvap(k)*otemp*oRv - 1.) & - + (-2.*lvap(k)*otemp*otemp*otemp*oRv) & - + otemp*otemp) - gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p - alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) & - * rvs_pp/rvs_p * rvs/rvs_p - alphsc = MAX(1.E-9, alphsc) - xsat = ssatw(k) - if (abs(xsat).lt. 1.E-9) xsat=0. - t1_evap = 2.*PI*( 1.0 - alphsc*xsat & - + 2.*alphsc*alphsc*xsat*xsat & - - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) & - / (1.+gamsc) - - Dc_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) & - * 4.*diffu(k)*ssatw(k)*rvs/rho_w) - idx_d = MAX(1, MIN(INT(1.E6*Dc_star), nbc)) - - idx_n = NINT(1.0 + FLOAT(nbc) * DLOG(nc(k)/t_Nc(1)) / nic1) - idx_n = MAX(1, MIN(idx_n, nbc)) - -!..Cloud water lookup table index. - if (rc(k).gt. r_c(1)) then - nic = NINT(ALOG10(rc(k))) - do nn = nic-1, nic+1 - n = nn - if ( (rc(k)/10.**nn).ge.1.0 .and. & - (rc(k)/10.**nn).lt.10.0) goto 159 - enddo - 159 continue - idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2) - idx_c = MAX(1, MIN(idx_c, ntb_c)) - else - idx_c = 1 - endif - - !prw_vcd(k) = MAX(DBLE(-rc(k)*orho*odt), & - ! -tpc_wev(idx_d, idx_c, idx_n)*orho*odt) - prw_vcd(k) = MAX(DBLE(-rc(k)*0.99*orho*odt), prw_vcd(k)) - pnc_wcd(k) = MAX(DBLE(-nc(k)*0.99*orho*odt), & - -tnc_wev(idx_d, idx_c, idx_n)*orho*odt) - endif else prw_vcd(k) = -rc(k)*orho*odt @@ -3104,7 +2622,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-IFDRY) rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k)) nc(k) = MAX(2., (nc1d(k) + DT*ncten(k))*rho(k)) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + nc(k) = Nt_c qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k)) temp(k) = t1d(k) + DT*tten(k) rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) @@ -3192,16 +2710,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) endif enddo -! #if ( WRF_CHEM == 1 ) -! do k = kts, kte -! evapprod(k) = prv_rev(k) - (min(zeroD0,prs_sde(k)) + & -! min(zeroD0,prg_gde(k))) -! rainprod(k) = prr_wau(k) + prr_rcw(k) + prs_scw(k) + & -! prg_scw(k) + prs_iau(k) + & -! prg_gcw(k) + prs_sci(k) + & -! pri_rci(k) -! enddo -! #endif !+---+-----------------------------------------------------------------+ !..Find max terminal fallspeed (distribution mass-weighted mean @@ -3278,10 +2786,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif enddo -!+---+-----------------------------------------------------------------+ - - if (.not. iiwarm) then - nstep = 0 do k = kte, kts, -1 vti = 0. @@ -3372,7 +2876,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & enddo if (ksed1(4) .eq. kte) ksed1(4) = kte-1 if (nstep .gt. 0) onstep(4) = 1./REAL(nstep) - endif !+---+-----------------------------------------------------------------+ !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD, @@ -3515,7 +3018,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !.. Instantly melt any cloud ice into cloud water if above 0C and !.. instantly freeze any cloud water found below HGFR. !+---+-----------------------------------------------------------------+ - if (.not. iiwarm) then do k = kts, kte xri = MAX(0.0, qi1d(k) + qiten(k)*DT) if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then @@ -3537,7 +3039,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY) endif enddo - endif !+---+-----------------------------------------------------------------+ !.. All tendencies computed, apply and pass back final values to parent. @@ -3609,8 +3110,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & end subroutine mp_thompson !+---+-----------------------------------------------------------------+ -!ctrlL -!+---+-----------------------------------------------------------------+ !..Creation of the lookup tables and support functions found below here. !+---+-----------------------------------------------------------------+ !..Rain collecting graupel (and inverse). Explicit CE integration. @@ -3620,24 +3119,16 @@ subroutine qr_acr_qg implicit none -!..Local variables INTEGER:: i, j, k, m, n, n2 INTEGER:: km, km_s, km_e DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2 - LOGICAL force_read_thompson, write_thompson_tables LOGICAL lexist,lopen INTEGER, allocatable :: good[:] - ! LOGICAL, EXTERNAL :: wrf_dm_on_monitor allocate(good[*]) -!+---+ - - ! CALL nl_get_force_read_thompson(1,force_read_thompson) - ! CALL nl_get_write_thompson_tables(1,write_thompson_tables) - good = 0 if (this_image()==1) then INQUIRE(FILE="qr_acr_qg.dat",EXIST=lexist) @@ -3660,12 +3151,6 @@ subroutine qr_acr_qg ! broadcast the data to all images do i=2,num_images() good[i] = good - ! tcg_racg(:,:,:,:)[i] = tcg_racg(:,:,:,:) - ! tmr_racg(:,:,:,:)[i] = tmr_racg(:,:,:,:) - ! tcr_gacr(:,:,:,:)[i] = tcr_gacr(:,:,:,:) - ! tmg_gacr(:,:,:,:)[i] = tmg_gacr(:,:,:,:) - ! tnr_racg(:,:,:,:)[i] = tnr_racg(:,:,:,:) - ! tnr_gacr(:,:,:,:)[i] = tnr_gacr(:,:,:,:) enddo ENDIF endif @@ -3695,12 +3180,8 @@ subroutine qr_acr_qg !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for !.. fortran indices. J. Michalakes, 2009Oct30. -! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) -! CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e ) -! #else km_s = 0 km_e = ntb_r*ntb_r1 - 1 -! #endif do km = km_s, km_e m = km / ntb_r1 + 1 @@ -3750,7 +3231,6 @@ subroutine qr_acr_qg z2 = z2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & *dvr*massg * N_g(n)* N_r(n2) enddo - 97 continue enddo tcg_racg(i,j,k,m) = t1 tmr_racg(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0) @@ -3764,15 +3244,6 @@ subroutine qr_acr_qg !..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30). -! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) -! CALL wrf_dm_gatherv(tcg_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tmr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tcr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tmg_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tnr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tnr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) -! #endif - IF ( this_image()==1 ) THEN print *, "Writing qr_acr_qg.dat in Thompson MP init" OPEN(63,file="qr_acr_qg.dat",form="unformatted",err=9234) @@ -3791,8 +3262,6 @@ subroutine qr_acr_qg end subroutine qr_acr_qg !+---+-----------------------------------------------------------------+ -!ctrlL -!+---+-----------------------------------------------------------------+ !..Rain collecting snow (and inverse). Explicit CE integration. !+---+-----------------------------------------------------------------+ @@ -3800,7 +3269,6 @@ subroutine qr_acr_qs implicit none -!..Local variables INTEGER:: i, j, k, m, n, n2 INTEGER:: km, km_s, km_e DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r @@ -3810,17 +3278,10 @@ subroutine qr_acr_qs DOUBLE PRECISION:: dvs, dvr, masss, massr DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4 DOUBLE PRECISION:: y1, y2, y3, y4 - LOGICAL force_read_thompson, write_thompson_tables LOGICAL lexist,lopen INTEGER, allocatable :: good[:] - ! LOGICAL, EXTERNAL :: wrf_dm_on_monitor allocate(good[*]) -!+---+ - - ! CALL nl_get_force_read_thompson(1,force_read_thompson) - ! CALL nl_get_write_thompson_tables(1,write_thompson_tables) - good = 0 IF ( this_image() == 1 ) THEN INQUIRE(FILE="qr_acr_qs.dat",EXIST=lexist) @@ -3847,18 +3308,6 @@ subroutine qr_acr_qs ENDIF do i=2,num_images() good[i] = good - ! tcs_racs1(:,:,:,:)[i] = tcs_racs1(:,:,:,:) - ! tmr_racs1(:,:,:,:)[i] = tmr_racs1(:,:,:,:) - ! tcs_racs2(:,:,:,:)[i] = tcs_racs2(:,:,:,:) - ! tmr_racs2(:,:,:,:)[i] = tmr_racs2(:,:,:,:) - ! tcr_sacr1(:,:,:,:)[i] = tcr_sacr1(:,:,:,:) - ! tms_sacr1(:,:,:,:)[i] = tms_sacr1(:,:,:,:) - ! tcr_sacr2(:,:,:,:)[i] = tcr_sacr2(:,:,:,:) - ! tms_sacr2(:,:,:,:)[i] = tms_sacr2(:,:,:,:) - ! tnr_racs1(:,:,:,:)[i] = tnr_racs1(:,:,:,:) - ! tnr_racs2(:,:,:,:)[i] = tnr_racs2(:,:,:,:) - ! tnr_sacr1(:,:,:,:)[i] = tnr_sacr1(:,:,:,:) - ! tnr_sacr2(:,:,:,:)[i] = tnr_sacr2(:,:,:,:) enddo ENDIF endif @@ -3882,7 +3331,6 @@ subroutine qr_acr_qs IF ( good .NE. 1 ) THEN if (this_image()==1) print *, "ThompMP: computing qr_acr_qs" do n2 = 1, nbr -! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) & + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) & - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2) @@ -3895,12 +3343,8 @@ subroutine qr_acr_qs !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for !.. fortran indices. J. Michalakes, 2009Oct30. -! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) -! CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e ) -! #else km_s = 0 km_e = ntb_r*ntb_r1 - 1 -! #endif do km = km_s, km_e m = km / ntb_r1 + 1 @@ -4033,21 +3477,6 @@ subroutine qr_acr_qs !..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30). -! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) -! CALL wrf_dm_gatherv(tcs_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tmr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tcs_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tmr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tcr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tms_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tcr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tms_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tnr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tnr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tnr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! CALL wrf_dm_gatherv(tnr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) -! #endif - IF ( this_image()==1 ) THEN print *, "Writing qr_acr_qs.dat in Thompson MP init" OPEN(63,file="qr_acr_qs.dat",form="unformatted",err=9234) @@ -4072,8 +3501,6 @@ subroutine qr_acr_qs end subroutine qr_acr_qs !+---+-----------------------------------------------------------------+ -!ctrlL -!+---+-----------------------------------------------------------------+ !..This is a literal adaptation of Bigg (1954) probability of drops of !..a particular volume freezing. Given this probability, simply freeze !..the proportion of drops summing their masses. @@ -4083,25 +3510,18 @@ subroutine freezeH2O implicit none -!..Local variables INTEGER:: i, j, k, m, n, n2 DOUBLE PRECISION, DIMENSION(nbr):: N_r, massr DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, & prob, vol, Texp, orho_w, & - lam_exp, lamr, N0_r, lamc, N0_c, y + lam_exp, lamr, N0_r, lamc, N0_c INTEGER:: nu_c REAL:: T_adjust - LOGICAL force_read_thompson, write_thompson_tables LOGICAL lexist,lopen INTEGER, allocatable :: good[:] - ! LOGICAL, EXTERNAL :: wrf_dm_on_monitor allocate(good[*]) -!+---+ - ! CALL nl_get_force_read_thompson(1,force_read_thompson) - ! CALL nl_get_write_thompson_tables(1,write_thompson_tables) - good = 0 IF ( this_image() == 1 ) THEN INQUIRE(FILE="freezeH2O.dat",EXIST=lexist) @@ -4122,12 +3542,6 @@ subroutine freezeH2O endif do i=2,num_images() good[i] = good - ! tpi_qrfz(:,:,:,:)[i] = tpi_qrfz(:,:,:,:) - ! tni_qrfz(:,:,:,:)[i] = tni_qrfz(:,:,:,:) - ! tpg_qrfz(:,:,:,:)[i] = tpg_qrfz(:,:,:,:) - ! tnr_qrfz(:,:,:,:)[i] = tnr_qrfz(:,:,:,:) - ! tpi_qcfz(:,:,:,:)[i] = tpi_qcfz(:,:,:,:) - ! tni_qcfz(:,:,:,:)[i] = tni_qcfz(:,:,:,:) enddo ENDIF ENDIF @@ -4159,7 +3573,6 @@ subroutine freezeH2O do m = 1, ntb_IN T_adjust = MAX(-3.0, MIN(3.0 - ALOG10(Nt_IN(m)), 3.0)) do k = 1, 45 -! print*, ' Freezing water for temp = ', -k Texp = DEXP( DFLOAT(k) - T_adjust*1.0D0 ) - 1.0D0 do j = 1, ntb_r1 do i = 1, ntb_r @@ -4230,8 +3643,6 @@ subroutine freezeH2O end subroutine freezeH2O !+---+-----------------------------------------------------------------+ -!ctrlL -!+---+-----------------------------------------------------------------+ !..Cloud ice converting to snow since portion greater than min snow !.. size. Given cloud ice content (kg/m**3), number concentration !.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into @@ -4246,14 +3657,11 @@ subroutine qi_aut_qs implicit none -!..Local variables INTEGER:: i, j, n2 DOUBLE PRECISION, DIMENSION(nbi):: N_i DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2 REAL:: xlimit_intg -!+---+ - do j = 1, ntb_i1 do i = 1, ntb_i lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi @@ -4286,7 +3694,6 @@ subroutine qi_aut_qs enddo end subroutine qi_aut_qs -!ctrlL !+---+-----------------------------------------------------------------+ !..Variable collision efficiency for rain collecting cloud water using !.. method of Beard and Grover, 1974 if a/A less than 0.25; otherwise @@ -4297,7 +3704,6 @@ subroutine table_Efrw implicit none -!..Local variables DOUBLE PRECISION:: vtr, stokes, reynolds, Ef_rw DOUBLE PRECISION:: p, yc0, F, G, H, z, K0, X INTEGER:: i, j @@ -4350,7 +3756,6 @@ subroutine table_Efrw enddo end subroutine table_Efrw -!ctrlL !+---+-----------------------------------------------------------------+ !..Variable collision efficiency for snow collecting cloud water using !.. method of Wang and Ji, 2000 except equate melted snow diameter to @@ -4361,7 +3766,6 @@ subroutine table_Efsw implicit none -!..Local variables DOUBLE PRECISION:: Ds_m, vts, vtc, stokes, reynolds, Ef_sw DOUBLE PRECISION:: p, yc0, F, G, H, z, K0 INTEGER:: i, j @@ -4394,7 +3798,6 @@ subroutine table_Efsw enddo end subroutine table_Efsw -!ctrlL !+---+-----------------------------------------------------------------+ !..Function to compute collision efficiency of collector species (rain, !.. snow, graupel) of aerosols. Follows Wang et al, 2010, ACP, which @@ -4439,7 +3842,6 @@ real function Eff_aero(D, Da, visc,rhoa,Temp,species) end function Eff_aero -!ctrlL !+---+-----------------------------------------------------------------+ !..Integrate rain size distribution from zero to D-star to compute the !.. number of drops smaller than D-star that evaporate in a single @@ -4451,13 +3853,10 @@ subroutine table_dropEvap implicit none -!..Local variables INTEGER:: i, j, k, n DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc DOUBLE PRECISION:: summ, summ2, lamc, N0_c INTEGER:: nu_c -! DOUBLE PRECISION:: Nt_r, N0, lam_exp, lam -! REAL:: xlimit_intg do n = 1, nbc massc(n) = am_r*Dc(n)**bm_r @@ -4469,16 +3868,13 @@ subroutine table_dropEvap lamc = (t_Nc(k)*am_r* ccg(2,nu_c)*ocg1(nu_c) / r_c(j))**obmr N0_c = t_Nc(k)*ocg1(nu_c) * lamc**cce(1,nu_c) do i = 1, nbc -!-GT tnc_wev(i,j,k) = GAMMP(nu_c+1., SNGL(Dc(i)*lamc))*t_Nc(k) N_c(i) = N0_c* Dc(i)**nu_c*EXP(-lamc*Dc(i))*dtc(i) -! if(j.eq.18 .and. k.eq.50) print*, ' N_c = ', N_c(i) summ = 0. summ2 = 0. do n = 1, i summ = summ + massc(n)*N_c(n) summ2 = summ2 + N_c(n) enddo -! if(j.eq.18 .and. k.eq.50) print*, ' DEBUG-TABLE: ', r_c(j), t_Nc(k), summ2, summ tpc_wev(i,j,k) = summ tnc_wev(i,j,k) = summ2 enddo @@ -4488,122 +3884,11 @@ subroutine table_dropEvap ! !..To do the same thing for rain. ! -! do k = 1, ntb_r -! do j = 1, ntb_r1 -! lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1 -! lam = lam_exp * (crg(3)*org2*org1)**obmr -! N0 = N0r_exp(j)/(crg(2)*lam_exp) * lam**cre(2) -! Nt_r = N0 * crg(2) / lam**cre(2) -! do i = 1, nbr -! xlimit_intg = lam*Dr(i) -! tnr_rev(i,j,k) = GAMMP(mu_r+1.0, xlimit_intg) * Nt_r -! enddo -! enddo -! enddo ! TO APPLY TABLE ABOVE !..Rain lookup table indexes. -! Dr_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) & -! * 0.78*4.*diffu(k)*xsat*rvs/rho_w) -! idx_d = NINT(1.0 + FLOAT(nbr) * DLOG(Dr_star/D0r) & -! / DLOG(Dr(nbr)/D0r)) -! idx_d = MAX(1, MIN(idx_d, nbr)) -! -! nir = NINT(ALOG10(rr(k))) -! do nn = nir-1, nir+1 -! n = nn -! if ( (rr(k)/10.**nn).ge.1.0 .and. & -! (rr(k)/10.**nn).lt.10.0) goto 154 -! enddo -!154 continue -! idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2) -! idx_r = MAX(1, MIN(idx_r, ntb_r)) -! -! lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr -! lam_exp = lamr * (crg(3)*org2*org1)**bm_r -! N0_exp = org1*rr(k)/am_r * lam_exp**cre(1) -! nir = NINT(DLOG10(N0_exp)) -! do nn = nir-1, nir+1 -! n = nn -! if ( (N0_exp/10.**nn).ge.1.0 .and. & -! (N0_exp/10.**nn).lt.10.0) goto 155 -! enddo -!155 continue -! idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3) -! idx_r1 = MAX(1, MIN(idx_r1, ntb_r1)) -! -! pnr_rev(k) = MIN(nr(k)*odts, SNGL(tnr_rev(idx_d,idx_r1,idx_r) & ! RAIN2M -! * odts)) end subroutine table_dropEvap -! -!ctrlL -!+---+-----------------------------------------------------------------+ -!..Fill the table of CCN activation data created from parcel model run -!.. by Trude Eidhammer with inputs of aerosol number concentration, -!.. vertical velocity, temperature, lognormal mean aerosol radius, and -!.. hygroscopicity, kappa. The data are read from external file and -!.. contain activated fraction of CCN for given conditions. -!+---+-----------------------------------------------------------------+ - - subroutine table_ccnAct - - ! USE module_domain - ! USE module_dm - implicit none - - ! LOGICAL, EXTERNAL:: wrf_dm_on_monitor - -!..Local variables - INTEGER:: iunit_mp_th1, i - LOGICAL:: opened - CHARACTER*64 errmess - - iunit_mp_th1 = -1 - ! IF ( wrf_dm_on_monitor() ) THEN - ! DO i = 20,99 - ! INQUIRE ( i , OPENED = opened ) - ! IF ( .NOT. opened ) THEN - ! iunit_mp_th1 = i - ! GOTO 2010 - ! ENDIF - ! ENDDO - ! 2010 CONTINUE - ! ENDIF -! #if defined(DM_PARALLEL) && !defined(STUBMPI) -! CALL wrf_dm_bcast_bytes ( iunit_mp_th1 , IWORDSIZE ) -! #endif - ! IF ( iunit_mp_th1 < 0 ) THEN - ! CALL wrf_error_fatal ( 'module_mp_thompson: table_ccnAct: '// & - ! 'Can not find unused fortran unit to read in lookup table.') - ! ENDIF - - ! IF ( wrf_dm_on_monitor() ) THEN - ! WRITE(errmess, '(A,I2)') 'module_mp_thompson: opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1 - ! CALL wrf_debug(150, errmess) - ! OPEN(iunit_mp_th1,FILE='CCN_ACTIVATE.BIN', & - ! FORM='UNFORMATTED',STATUS='OLD',ERR=9009) - ! ENDIF - -! #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes(A, size(A)*R4SIZE) - -! IF ( wrf_dm_on_monitor() ) READ(iunit_mp_th1,ERR=9010) tnccn_act -! #if defined(DM_PARALLEL) && !defined(STUBMPI) -! DM_BCAST_MACRO(tnccn_act) -! #endif - - - RETURN - 9009 CONTINUE - ! WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1 - ! CALL wrf_error_fatal(errmess) - RETURN - 9010 CONTINUE - ! WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error reading CCN_ACTIVATE.BIN on unit ',iunit_mp_th1 - ! CALL wrf_error_fatal(errmess) - - end subroutine table_ccnAct -!^L !+---+-----------------------------------------------------------------+ !..Retrieve fraction of CCN that gets activated given the model temp, !.. vertical velocity, and available CCN concentration. The lookup @@ -4622,13 +3907,6 @@ real function activ_ncloud(Tt, Ww, NCCN) INTEGER:: i, j, k, l, m, n REAL:: A, B, C, D, t, u, x1, x2, y1, y2, nx, wy, fraction - -! ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) ntb_arc -! ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) ntb_arw -! ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) ntb_art -! ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) ntb_arr -! ta_Ka = (/0.2, 0.4, 0.6, 0.8/) ntb_ark - n_local = NCCN * 1.E-6 w_local = Ww @@ -4678,21 +3956,12 @@ real function activ_ncloud(Tt, Ww, NCCN) t = (nx-x1)/(x2-x1) u = (wy-y1)/(y2-y1) -! t = (n_local-ta(Na(i-1))/(ta_Na(i)-ta_Na(i-1)) -! u = (w_local-ta_Ww(j-1))/(ta_Ww(j)-ta_Ww(j-1)) - fraction = (1.0-t)*(1.0-u)*A + t*(1.0-u)*B + t*u*C + (1.0-t)*u*D -! if (NCCN*fraction .gt. 0.75*Nt_c_max) then -! write(*,*) ' DEBUG-GT ', n_local, w_local, Tt, i, j, k -! endif - activ_ncloud = NCCN*fraction end function activ_ncloud -!+---+-----------------------------------------------------------------+ -!+---+-----------------------------------------------------------------+ SUBROUTINE GCF(GAMMCF,A,X,GLN) ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS ! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS @@ -4836,7 +4105,6 @@ REAL FUNCTION RSLF(P,T) X=MAX(-80.,T-273.16) -! ESL=612.2*EXP(17.67*X/(T-29.65)) ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) ESL=MIN(ESL, P*0.15) ! Even with P=1050mb and T=55C, the sat. vap. pres only contributes to ~15% of total pres. RSLF=.622*ESL/(P-ESL) @@ -4882,71 +4150,19 @@ REAL FUNCTION RSIF(P,T) END FUNCTION RSIF -!+---+-----------------------------------------------------------------+ - real function iceDeMott(tempc, qv, qvs, qvsi, rho, nifa) + real function iceDeMott(tempc, rho, nifa) implicit none - REAL, INTENT(IN):: tempc, qv, qvs, qvsi, rho, nifa - -!..Local vars - REAL:: satw, sati, siw, p_x, si0x, dtt, dsi, dsw, dab, fc, hx - REAL:: ntilde, n_in, nmax, nhat, mux, xni, nifa_cc - REAL, PARAMETER:: p_c1 = 1000. - REAL, PARAMETER:: p_rho_c = 0.76 - REAL, PARAMETER:: p_alpha = 1.0 - REAL, PARAMETER:: p_gam = 2. - REAL, PARAMETER:: delT = 5. - REAL, PARAMETER:: T0x = -40. - REAL, PARAMETER:: Sw0x = 0.97 - REAL, PARAMETER:: delSi = 0.1 - REAL, PARAMETER:: hdm = 0.15 - REAL, PARAMETER:: p_psi = 0.058707*p_gam/p_rho_c - REAL, PARAMETER:: aap = 1. - REAL, PARAMETER:: bbp = 0. - REAL, PARAMETER:: y1p = -35. - REAL, PARAMETER:: y2p = -25. - REAL, PARAMETER:: rho_not0 = 101325./(287.05*273.15) + REAL, INTENT(IN):: tempc, rho, nifa -!+---+ + REAL:: xni, nifa_cc + REAL, PARAMETER:: rho_not0 = 101325./(287.05*273.15) xni = 0.0 -! satw = qv/qvs -! sati = qv/qvsi -! siw = qvs/qvsi -! p_x = -1.0261+(3.1656e-3*tempc)+(5.3938e-4*(tempc*tempc)) & -! + (8.2584e-6*(tempc*tempc*tempc)) -! si0x = 1.+(10.**p_x) -! if (sati.ge.si0x .and. satw.lt.0.985) then -! dtt = delta_p (tempc, T0x, T0x+delT, 1., hdm) -! dsi = delta_p (sati, Si0x, Si0x+delSi, 0., 1.) -! dsw = delta_p (satw, Sw0x, 1., 0., 1.) -! fc = dtt*dsi*0.5 -! hx = min(fc+((1.-fc)*dsw), 1.) -! ntilde = p_c1*p_gam*((exp(12.96*(sati-1.1)))**0.3) / p_rho_c -! if (tempc .le. y1p) then -! n_in = ntilde -! elseif (tempc .ge. y2p) then -! n_in = p_psi*p_c1*exp(12.96*(sati-1.)-0.639) -! else -! if (tempc .le. -30.) then -! nmax = p_c1*p_gam*(exp(12.96*(siw-1.1)))**0.3/p_rho_c -! else -! nmax = p_psi*p_c1*exp(12.96*(siw-1.)-0.639) -! endif -! ntilde = MIN(ntilde, nmax) -! nhat = MIN(p_psi*p_c1*exp(12.96*(sati-1.)-0.639), nmax) -! dab = delta_p (tempc, y1p, y2p, aap, bbp) -! n_in = MIN(nhat*(ntilde/nhat)**dab, nmax) -! endif -! mux = hx*p_alpha*n_in*rho -! xni = mux*((6700.*nifa)-200.)/((6700.*5.E5)-200.) -! elseif (satw.ge.0.985 .and. tempc.gt.HGFR-273.15) then nifa_cc = nifa*RHO_NOT0*1.E-6/rho -! xni = 3.*nifa_cc**(1.25)*exp((0.46*(-tempc))-11.6) ! [DeMott, 2015] xni = (5.94e-5*(-tempc)**3.33) & ! [DeMott, 2010] * (nifa_cc**((-0.0264*(tempc))+0.0033)) xni = xni*rho/RHO_NOT0 * 1000. -! endif iceDeMott = MAX(0., xni) @@ -5018,9 +4234,6 @@ REAL FUNCTION delta_p (yy, y1, y2, aa, bb) END FUNCTION delta_p -!+---+-----------------------------------------------------------------+ -!ctrlL - !+---+-----------------------------------------------------------------+ !..Compute _radiation_ effective radii of cloud water, ice, and snow. !.. These are entirely consistent with microphysics assumptions, not @@ -5035,12 +4248,11 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & IMPLICIT NONE -!..Sub arguments INTEGER, INTENT(IN):: kts, kte REAL, DIMENSION(kts:kte), INTENT(IN):: & & t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: re_qc1d, re_qi1d, re_qs1d -!..Local variables + INTEGER:: k REAL, DIMENSION(kts:kte):: rho, rc, nc, ri, ni, rs REAL:: smo2, smob, smoc @@ -5059,7 +4271,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) rc(k) = MAX(R1, qc1d(k)*rho(k)) nc(k) = MAX(R2, nc1d(k)*rho(k)) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + nc(k) = Nt_c if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true. ri(k) = MAX(R1, qi1d(k)*rho(k)) ni(k) = MAX(R2, ni1d(k)*rho(k)) @@ -5133,319 +4345,4 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & end subroutine calc_effectRad -!+---+-----------------------------------------------------------------+ -!..Compute radar reflectivity assuming 10 cm wavelength radar and using -!.. Rayleigh approximation. Only complication is melted snow/graupel -!.. which we treat as water-coated ice spheres and use Uli Blahak's -!.. library of routines. The meltwater fraction is simply the amount -!.. of frozen species remaining from what initially existed at the -!.. melting level interface. -!+---+-----------------------------------------------------------------+ - -! subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & -! t1d, p1d, dBZ, kts, kte, ii, jj) -! -! IMPLICIT NONE -! -! !..Sub arguments -! INTEGER, INTENT(IN):: kts, kte, ii, jj -! REAL, DIMENSION(kts:kte), INTENT(IN):: & -! qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d -! REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ -! ! REAL, DIMENSION(kts:kte), INTENT(INOUT):: vt_dBZ -! -! !..Local variables -! REAL, DIMENSION(kts:kte):: temp, pres, qv, rho, rhof -! REAL, DIMENSION(kts:kte):: rc, rr, nr, rs, rg -! -! DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g -! REAL, DIMENSION(kts:kte):: mvd_r -! REAL, DIMENSION(kts:kte):: smob, smo2, smoc, smoz -! REAL:: oM3, M0, Mrat, slam1, slam2, xDs -! REAL:: ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts -! REAL:: vtr_dbz_wt, vts_dbz_wt, vtg_dbz_wt -! -! REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel -! -! DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg -! REAL:: a_, b_, loga_, tc0 -! DOUBLE PRECISION:: fmelt_s, fmelt_g -! -! INTEGER:: i, k, k_0, kbot, n -! LOGICAL:: melti -! LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg -! -! DOUBLE PRECISION:: cback, x, eta, f_d -! REAL:: xslw1, ygra1, zans1 -! -! !+---+ -! -! do k = kts, kte -! dBZ(k) = -35.0 -! enddo -! -! !+---+-----------------------------------------------------------------+ -! !..Put column of data into local arrays. -! !+---+-----------------------------------------------------------------+ -! do k = kts, kte -! temp(k) = t1d(k) -! qv(k) = MAX(1.E-10, qv1d(k)) -! pres(k) = p1d(k) -! rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) -! rhof(k) = SQRT(RHO_NOT/rho(k)) -! rc(k) = MAX(R1, qc1d(k)*rho(k)) -! if (qr1d(k) .gt. R1) then -! rr(k) = qr1d(k)*rho(k) -! nr(k) = MAX(R2, nr1d(k)*rho(k)) -! lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr -! ilamr(k) = 1./lamr -! N0_r(k) = nr(k)*org2*lamr**cre(2) -! mvd_r(k) = (3.0 + mu_r + 0.672) * ilamr(k) -! L_qr(k) = .true. -! else -! rr(k) = R1 -! nr(k) = R1 -! mvd_r(k) = 50.E-6 -! L_qr(k) = .false. -! endif -! if (qs1d(k) .gt. R2) then -! rs(k) = qs1d(k)*rho(k) -! L_qs(k) = .true. -! else -! rs(k) = R1 -! L_qs(k) = .false. -! endif -! if (qg1d(k) .gt. R2) then -! rg(k) = qg1d(k)*rho(k) -! L_qg(k) = .true. -! else -! rg(k) = R1 -! L_qg(k) = .false. -! endif -! enddo -! -! !+---+-----------------------------------------------------------------+ -! !..Calculate y-intercept, slope, and useful moments for snow. -! !+---+-----------------------------------------------------------------+ -! do k = kts, kte -! tc0 = MIN(-0.1, temp(k)-273.15) -! smob(k) = rs(k)*oams -! -! !..All other moments based on reference, 2nd moment. If bm_s.ne.2, -! !.. then we must compute actual 2nd moment and use as reference. -! if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then -! smo2(k) = smob(k) -! else -! loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & -! & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & -! & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & -! & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & -! & + sa(10)*bm_s*bm_s*bm_s -! a_ = 10.0**loga_ -! b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & -! & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & -! & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & -! & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & -! & + sb(10)*bm_s*bm_s*bm_s -! smo2(k) = (smob(k)/a_)**(1./b_) -! endif -! -! !..Calculate bm_s+1 (th) moment. Useful for diameter calcs. -! loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & -! & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & -! & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & -! & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & -! & + sa(10)*cse(1)*cse(1)*cse(1) -! a_ = 10.0**loga_ -! b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & -! & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & -! & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & -! & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) -! smoc(k) = a_ * smo2(k)**b_ -! -! !..Calculate bm_s*2 (th) moment. Useful for reflectivity. -! loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(3) & -! & + sa(4)*tc0*cse(3) + sa(5)*tc0*tc0 & -! & + sa(6)*cse(3)*cse(3) + sa(7)*tc0*tc0*cse(3) & -! & + sa(8)*tc0*cse(3)*cse(3) + sa(9)*tc0*tc0*tc0 & -! & + sa(10)*cse(3)*cse(3)*cse(3) -! a_ = 10.0**loga_ -! b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(3) + sb(4)*tc0*cse(3) & -! & + sb(5)*tc0*tc0 + sb(6)*cse(3)*cse(3) & -! & + sb(7)*tc0*tc0*cse(3) + sb(8)*tc0*cse(3)*cse(3) & -! & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(3)*cse(3)*cse(3) -! smoz(k) = a_ * smo2(k)**b_ -! enddo -! -! !+---+-----------------------------------------------------------------+ -! !..Calculate y-intercept, slope values for graupel. -! !+---+-----------------------------------------------------------------+ -! -! N0_min = gonv_max -! k_0 = kts -! do k = kte, kts, -1 -! if (temp(k).ge.270.65) k_0 = MAX(k_0, k) -! enddo -! do k = kte, kts, -1 -! if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then -! xslw1 = 4.01 + alog10(mvd_r(k)) -! else -! xslw1 = 0.01 -! endif -! ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) -! zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1)) -! N0_exp = 10.**(zans1) -! N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) -! N0_min = MIN(N0_exp, N0_min) -! N0_exp = N0_min -! lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 -! lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg -! ilamg(k) = 1./lamg -! N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) -! enddo -! -! !+---+-----------------------------------------------------------------+ -! !..Locate K-level of start of melting (k_0 is level above). -! !+---+-----------------------------------------------------------------+ -! melti = .false. -! k_0 = kts -! do k = kte-1, kts, -1 -! if ( (temp(k).gt.273.15) .and. L_qr(k) & -! & .and. (L_qs(k+1).or.L_qg(k+1)) ) then -! k_0 = MAX(k+1, k_0) -! melti=.true. -! goto 195 -! endif -! enddo -! 195 continue -! -! !+---+-----------------------------------------------------------------+ -! !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps) -! !.. and non-water-coated snow and graupel when below freezing are -! !.. simple. Integrations of m(D)*m(D)*N(D)*dD. -! !+---+-----------------------------------------------------------------+ -! -! do k = kts, kte -! ze_rain(k) = 1.e-22 -! ze_snow(k) = 1.e-22 -! ze_graupel(k) = 1.e-22 -! if (L_qr(k)) ze_rain(k) = N0_r(k)*crg(4)*ilamr(k)**cre(4) -! if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & -! & * (am_s/900.0)*(am_s/900.0)*smoz(k) -! if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & -! & * (am_g/900.0)*(am_g/900.0) & -! & * N0_g(k)*cgg(4)*ilamg(k)**cge(4) -! enddo -! -! !+---+-----------------------------------------------------------------+ -! !..Special case of melting ice (snow/graupel) particles. Assume the -! !.. ice is surrounded by the liquid water. Fraction of meltwater is -! !.. extremely simple based on amount found above the melting level. -! !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting -! !.. routines). -! !+---+-----------------------------------------------------------------+ -! -! if (.not. iiwarm .and. melti .and. k_0.ge.2) then -! do k = k_0-1, kts, -1 -! -! !..Reflectivity contributed by melting snow -! if (L_qs(k) .and. L_qs(k_0) ) then -! fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) -! eta = 0.d0 -! oM3 = 1./smoc(k) -! M0 = (smob(k)*oM3) -! Mrat = smob(k)*M0*M0*M0 -! slam1 = M0 * Lam0 -! slam2 = M0 * Lam1 -! do n = 1, nrbins -! x = am_s * xxDs(n)**bm_s -! call rayleigh_soak_wetgraupel (x, DBLE(ocms), DBLE(obms), & -! & fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, & -! & CBACK, mixingrulestring_s, matrixstring_s, & -! & inclusionstring_s, hoststring_s, & -! & hostmatrixstring_s, hostinclusionstring_s) -! f_d = Mrat*(Kap0*DEXP(-slam1*xxDs(n)) & -! & + Kap1*(M0*xxDs(n))**mu_s * DEXP(-slam2*xxDs(n))) -! eta = eta + f_d * CBACK * simpson(n) * xdts(n) -! enddo -! ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta) -! endif -! -! !..Reflectivity contributed by melting graupel -! -! if (L_qg(k) .and. L_qg(k_0) ) then -! fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) -! eta = 0.d0 -! lamg = 1./ilamg(k) -! do n = 1, nrbins -! x = am_g * xxDg(n)**bm_g -! call rayleigh_soak_wetgraupel (x, DBLE(ocmg), DBLE(obmg), & -! & fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, & -! & CBACK, mixingrulestring_g, matrixstring_g, & -! & inclusionstring_g, hoststring_g, & -! & hostmatrixstring_g, hostinclusionstring_g) -! f_d = N0_g(k)*xxDg(n)**mu_g * DEXP(-lamg*xxDg(n)) -! eta = eta + f_d * CBACK * simpson(n) * xdtg(n) -! enddo -! ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta) -! endif -! -! enddo -! endif -! -! do k = kte, kts, -1 -! dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18) -! enddo -! -! -! !..Reflectivity-weighted terminal velocity (snow, rain, graupel, mix). -! ! do k = kte, kts, -1 -! ! vt_dBZ(k) = 1.E-3 -! ! if (rs(k).gt.R2) then -! ! Mrat = smob(k) / smoc(k) -! ! ils1 = 1./(Mrat*Lam0 + fv_s) -! ! ils2 = 1./(Mrat*Lam1 + fv_s) -! ! t1_vts = Kap0*csg(5)*ils1**cse(5) -! ! t2_vts = Kap1*Mrat**mu_s*csg(11)*ils2**cse(11) -! ! ils1 = 1./(Mrat*Lam0) -! ! ils2 = 1./(Mrat*Lam1) -! ! t3_vts = Kap0*csg(6)*ils1**cse(6) -! ! t4_vts = Kap1*Mrat**mu_s*csg(12)*ils2**cse(12) -! ! vts_dbz_wt = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) -! ! if (temp(k).ge.273.15 .and. temp(k).lt.275.15) then -! ! vts_dbz_wt = vts_dbz_wt*1.5 -! ! elseif (temp(k).ge.275.15) then -! ! vts_dbz_wt = vts_dbz_wt*2.0 -! ! endif -! ! else -! ! vts_dbz_wt = 1.E-3 -! ! endif -! -! ! if (rr(k).gt.R1) then -! ! lamr = 1./ilamr(k) -! ! vtr_dbz_wt = rhof(k)*av_r*crg(13)*(lamr+fv_r)**(-cre(13)) & -! ! & / (crg(4)*lamr**(-cre(4))) -! ! else -! ! vtr_dbz_wt = 1.E-3 -! ! endif -! -! ! if (rg(k).gt.R2) then -! ! lamg = 1./ilamg(k) -! ! vtg_dbz_wt = rhof(k)*av_g*cgg(5)*lamg**(-cge(5)) & -! ! & / (cgg(4)*lamg**(-cge(4))) -! ! else -! ! vtg_dbz_wt = 1.E-3 -! ! endif -! -! ! vt_dBZ(k) = (vts_dbz_wt*ze_snow(k) + vtr_dbz_wt*ze_rain(k) & -! ! & + vtg_dbz_wt*ze_graupel(k)) & -! ! & / (ze_rain(k)+ze_snow(k)+ze_graupel(k)) -! ! enddo -! -! end subroutine calc_refl10cm -! -!+---+-----------------------------------------------------------------+ - -!+---+-----------------------------------------------------------------+ END MODULE module_mp_thompson -!+---+-----------------------------------------------------------------+ diff --git a/test/test-initialization.f90 b/test/test-initialization.f90 index 6da5b589..837f7b23 100644 --- a/test/test-initialization.f90 +++ b/test/test-initialization.f90 @@ -24,7 +24,7 @@ program main block type(domain_t) :: domain - integer :: i,nz, ypos + integer :: i,nz real :: start, finish call cpu_time(start)