#define MODAL_AERO MODULE module_bl_camuwpbl_driver !Note: Comments starting with "!!" are directly taken from CAM's interface routine for the UW PBL (vertical_diffusion.F90) !This modules is based on vertical_diffusion.F90 of CAM5 !List of possible modifications in the future: !_____________________________________________ !1. Constituents(per Kg of DRY air or "dry" constituents) are not diffused currently. The call which ! makes DRY constituents diffusion decision flag 'true' is commented out in camuwpblinit subroutine. ! Therefore the flag is ALWAYS false for DRY constituent diffusion. ! ! REASON: DRY constituents are Aerosols and other species (excepts for cloud mass mixing ratios and number concentrations). ! ------ In WRF, DRY constituents appear only when WRF_CHEM simulations are run. DRY constituents are vertically mixed(diffused) ! in dry_dep_driver of WRF_CHEM. Therefore, DRY constituents are not diffused in this PBL subroutine ! !2. The liquid number concentrations is NOT diffused(or mixed) in PBL to mimic CAM5, which diffuses liquid number ! concentrations in dropmixnuc. !3. Surface fluxes for ALL the constituents are ZERO except for the water vapours (1st constituent) !4. Molecular diffusion is turned off for now !5. Mountain stresses are not computed for now (all calls to Mountain stresses currently may hold undefined variables) !6. *Formulas used for computing surface stresses is based on the formula given in CAM5's BareGround.F90. ! This formula may not work well for ocean surfaces (possible future modification) !_____________________________________________ !!----------------------------------------------------------------------------------------------------- ! !! Module to compute vertical diffusion of momentum, moisture, trace constituents ! !! and static energy. Separate modules compute ! !! 1. stresses associated with turbulent flow over orography ! !! ( turbulent mountain stress ) ! !! 2. eddy diffusivities, including nonlocal tranport terms ! !! 3. molecular diffusivities ! !! 4. coming soon... gravity wave drag ! !! Lastly, a implicit diffusion solver is called, and tendencies retrieved by ! !! differencing the diffused and initial states. ! !! ! !! Calling sequence: ! !! ! !! vertical_diffusion_init Initializes vertical diffustion constants and modules ! !! init_molec_diff Initializes molecular diffusivity module ! !! init_eddy_diff Initializes eddy diffusivity module (includes PBL) ! !! init_tms Initializes turbulent mountain stress module ! !! init_vdiff Initializes diffusion solver module ! !! vertical_diffusion_ts_init Time step initialization (only used for upper boundary condition) ! !! vertical_diffusion_tend Computes vertical diffusion tendencies ! !! compute_tms Computes turbulent mountain stresses ! !! compute_eddy_diff Computes eddy diffusivities and countergradient terms ! !! compute_vdiff Solves vertical diffusion equations, including molecular diffusivities ! !! ! !!---------------------------Code history-------------------------------------------------------------- ! !! J. Rosinski : Jun. 1992 ! !! J. McCaa : Sep. 2004 ! !! S. Park : Aug. 2006, Dec. 2008. Jan. 2010 ! ! B. Singh : Nov. 2010 (ported to WRF by balwinder.singh@pnl.gov) !!----------------------------------------------------------------------------------------------------- ! use shr_kind_mod, only : r8 => shr_kind_r8 use module_cam_support, only : pcols, pver, pverp, endrun, iulog,fieldname_len,pcnst =>pcnst_runtime use constituents, only : qmin use diffusion_solver, only : vdiff_selector use physconst, only : & cpair , & ! Specific heat of dry air gravit , & ! Acceleration due to gravity rair , & ! Gas constant for dry air zvir , & ! rh2o/rair - 1 latvap , & ! Latent heat of vaporization latice , & ! Latent heat of fusion karman , & ! von Karman constant mwdry , & ! Molecular weight of dry air avogad , & ! Avogadro's number boltz ! Boltzman's constant implicit none private save !! ----------------- ! !! Public interfaces ! !! ----------------- ! public camuwpblinit ! Initialization public camuwpbl ! Driver for the PBL scheme public vd_register ! Init routine for constituents !! ------------ ! !! Private data ! !! ------------ ! character(len=16) :: eddy_scheme !! Default set in phys_control.F90, use namelist to change !! 'HB' = Holtslag and Boville (default) !! 'HBR' = Holtslag and Boville and Rash !! 'diag_TKE' = Bretherton and Park ( UW Moist Turbulence Scheme ) integer, parameter :: nturb = 5 !! Number of iterations for solution ( when 'diag_TKE' scheme is selected ) logical, parameter :: wstarent = .true. !! Use wstar (.true.) or TKE (.false.) entrainment closure ( when 'diag_TKE' scheme is selected ) logical :: do_pseudocon_diff = .false. !! If .true., do pseudo-conservative variables diffusion character(len=16) :: shallow_scheme !! For checking compatibility between eddy diffusion and shallow convection schemes !! 'Hack' = Hack Shallow Convection Scheme !! 'UW' = Park and Bretherton ( UW Shallow Convection Scheme ) character(len=16) :: microp_scheme !! Microphysics scheme logical :: do_molec_diff = .false. !! Switch for molecular diffusion logical :: do_tms !! Switch for turbulent mountain stress real(r8) :: tms_orocnst !! Converts from standard deviation to height real(r8) :: tms_z0fac ! Converts from standard deviation to height type(vdiff_selector) :: fieldlist_wet !! Logical switches for moist mixing ratio diffusion type(vdiff_selector) :: fieldlist_dry !! Logical switches for dry mixing ratio diffusion integer :: ntop !! Top interface level to which vertical diffusion is applied ( = 1 ). integer :: nbot !! Bottom interface level to which vertical diffusion is applied ( = pver ). integer :: tke_idx, kvh_idx, kvm_idx !! TKE and eddy diffusivity indices for fields in the physics buffer integer :: turbtype_idx, smaw_idx !! Turbulence type and instability functions integer :: tauresx_idx, tauresy_idx !! Redisual stress for implicit surface stress integer :: ixcldice, ixcldliq !! Constituent indices for cloud liquid and ice water integer :: ixnumice, ixnumliq #ifdef MODAL_AERO integer :: ixndrop #endif integer :: wgustd_index CONTAINS subroutine camuwpbl(dt,u_phy,v_phy,th_phy,rho,qv_curr,hfx,qfx,ustar,p8w & ,p_phy,z,t_phy,qc_curr,qi_curr,z_at_w,cldfra_old_mp,cldfra,ht & ,rthratenlw,exner,is_CAMMGMP_used & ,itimestep,qnc_curr,qni_curr,wsedl3d & ,ids,ide, jds,jde, kds,kde & ,ims,ime, jms,jme, kms,kme & ,its,ite, jts,jte, kts,kte & !Intent-inout ,tauresx2d,tauresy2d & ,rublten,rvblten,rthblten,rqiblten,rqniblten,rqvblten,rqcblten & ,kvm3d,kvh3d & !Intent-out ,tpert2d,qpert2d,wpert2d,smaw3d,turbtype3d & ,tke_pbl,pblh2d,kpbl2d ) !!---------------------------------------------------- ! !! This is an interface routine for vertical diffusion ! ! This subroutine is called by module_pbl_driver and ! ! it calls: wrf_error_fatal,compute_tms, ! ! compute_eddy_diff,aqsat,compute_vdiff and ! ! positive_moisture. !!---------------------------------------------------- ! use module_cam_support, only : pcols use trb_mtn_stress, only : compute_tms use eddy_diff, only : compute_eddy_diff use wv_saturation, only : fqsatd, aqsat use molec_diff, only : compute_molec_diff, vd_lu_qdecomp use constituents, only : qmincg, qmin use diffusion_solver !!, only : compute_vdiff, any, operator(.not.) #ifdef MODAL_AERO use modal_aero_data #endif implicit none !------------------------------------------------------------------------! ! Input ! !------------------------------------------------------------------------! logical, intent(in) :: is_CAMMGMP_used integer, intent(in) :: ids,ide, jds,jde, kds,kde integer, intent(in) :: ims,ime, jms,jme, kms,kme integer, intent(in) :: its,ite, jts,jte, kts,kte integer, intent(in) :: itimestep real, intent(in) :: dt ! Time step (s) real, dimension( ims:ime,jms:jme ), intent(in) :: hfx !Sensible heat flux at surface (w/m2) real, dimension( ims:ime,jms:jme ), intent(in) :: qfx !Moisture flux at surface (kg m-2 s-1) real, dimension( ims:ime,jms:jme ), intent(in) :: ustar !Friction velocity (m/s) real, dimension( ims:ime,jms:jme ), intent(in) :: ht !Terrain height (m) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: v_phy !Y-component of wind (m/s) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: u_phy !X-component of wind (m/s) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: th_phy !Potential temperature (K) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rho !Air density (kg/m3) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qv_curr !Water vapor mixing ratio - Moisture (kg/kg) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qc_curr !Cloud water mixing ratio - Cloud liq (kg/kg) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qi_curr !Ice mixing ratio - Cloud ice (kg/kg) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qnc_curr !Liq # mixing ratio - Cloud liq # (#/kg) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qni_curr !Ice # mixing ratio - Cloud ice # (#/kg) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p_phy !Pressure at mid-level (Pa) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p8w !Pressure at level interface (Pa) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z !Height above sea level at mid-level (m) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: t_phy !Temperature (K) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z_at_w !Height above sea level at layer interfaces (m) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cldfra_old_mp !Cloud fraction [unitless] real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cldfra !Cloud fraction [unitless] real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: exner !Dimensionless pressure [unitless] real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rthratenlw !Tendency for LW ( K/s) real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: wsedl3d !Sedimentation velocity of stratiform liquid cloud droplet (m/s) !------------------------------------------------------------------------! ! Output ! !------------------------------------------------------------------------! real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rublten !Tendency for u_phy (Pa m s-2) real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rvblten !Tendency for v_phy (Pa m s-2) real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rthblten !Tendency for th_phy (Pa K s-1) real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqvblten !Tendency for qv_curr (Pa kg kg-1 s-1) real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqcblten !Tendency for qc_curr (Pa kg kg-1 s-1) real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqiblten !Tendency for qi_curr (Pa kg kg-1 s-1) real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqniblten !Tendency for qni_curr(Pa # kg-1 s-1) integer, dimension( ims:ime,jms:jme ), intent(out) :: kpbl2d !Layer index containing PBL top within or at the base interface real, dimension( ims:ime,jms:jme ), intent(out) :: pblh2d !Planetary boundary layer height (m) real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: tke_pbl !Turbulence kinetic energy at midpoints (m^2/s^2) real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: turbtype3d !Turbulent interface types [ no unit ] real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: smaw3d !Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units] !---------------------------------------------------------------------------! ! Local, carried on from one timestep to the other (defined in registry)! !---------------------------------------------------------------------------! real, dimension( ims:ime, jms:jme ) , intent(inout ):: tauresx2d,tauresy2d !X AND Y-COMP OF RESIDUAL STRESSES(m^2/s^2) real, dimension( ims:ime, jms:jme ) , intent(out) :: tpert2d ! Convective temperature excess (K) real, dimension( ims:ime, jms:jme ) , intent(out) :: qpert2d ! Convective humidity excess (kg/kg) real, dimension( ims:ime, jms:jme ) , intent(out) :: wpert2d ! Turbulent velocity excess (m/s) real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: kvm3d,kvh3d !Eddy diffusivity for momentum and heat(m^2/s) !---------------------------------------------------------------------------! ! Local ! !---------------------------------------------------------------------------! character(128) :: errstring ! Error status for compute_vdiff logical :: kvinit ! Tell compute_eddy_diff/ caleddy to initialize kvh, kvm (uses kvf) logical :: is_first_step ! Flag to know if this a first time step integer :: i,j,k,itsp1,itile_len,ktep1,kflip,ncol,ips,icnst,m,kp1 integer :: lchnk ! Chunk identifier real(r8) :: tauFac, uMean, dp, multFrc real(r8) :: ztodt ! 2*delta-t (s) real(r8) :: rztodt ! 1./ztodt (1/s) real(r8),parameter :: invalidVal = -999888777.0_r8 real(r8) :: topflx( pcols) ! Molecular heat flux at top interface real(r8) :: wpert( pcols) ! Turbulent velocity excess (m/s) real(r8) :: tauresx( pcols) ! [Residual stress to be added in vdiff to correct... real(r8) :: tauresy( pcols) ! for turb stress mismatch between sfc and atm accumulated.] (N/m2) real(r8) :: ipbl( pcols) ! If 1, PBL is CL, while if 0, PBL is STL. real(r8) :: kpblh( pcols) ! Layer index containing PBL top within or at the base interface real(r8) :: wstarPBL(pcols) ! Convective velocity within PBL (m/s) real(r8) :: sgh( pcols) ! Standard deviation of orography (m) real(r8) :: landfrac(pcols) ! Land fraction [unitless] real(r8) :: taux( pcols) ! x surface stress (N/m2) real(r8) :: tauy( pcols) ! y surface stress (N/m2) real(r8) :: tautotx( pcols) ! U component of total surface stress (N/m2) real(r8) :: tautoty( pcols) ! V component of total surface stress (N/m2) real(r8) :: ksrftms( pcols) ! Turbulent mountain stress surface drag coefficient (kg/s/m2) real(r8) :: tautmsx( pcols) ! U component of turbulent mountain stress (N/m2) real(r8) :: tautmsy( pcols) ! V component of turbulent mountain stress (N/m2) real(r8) :: ustar8( pcols) ! Surface friction velocity (m/s) real(r8) :: pblh( pcols) ! Planetary boundary layer height (m) real(r8) :: tpert( pcols) ! Convective temperature excess (K) real(r8) :: qpert( pcols) ! Convective humidity excess (kg/kg) real(r8) :: shflx( pcols) ! Surface sensible heat flux (w/m2) real(r8) :: phis( pcols) ! Geopotential at terrain height (m2/s2) real(r8) :: cldn8( pcols,kte) ! New stratus fraction (fraction) real(r8) :: qrl8( pcols,kte) ! LW radiative cooling rate(W/kg*Pa) real(r8) :: wsedl8( pcols,kte) ! Sedimentation velocity of stratiform liquid cloud droplet (m/s) real(r8) :: dtk( pcols,kte) ! T tendency from KE dissipation real(r8) :: qt( pcols,kte) ! real(r8) :: sl_prePBL( pcols,kte) ! real(r8) :: qt_prePBL( pcols,kte) ! real(r8) :: slten( pcols,kte) ! real(r8) :: qtten( pcols,kte) ! real(r8) :: sl( pcols,kte) ! real(r8) :: ftem( pcols,kte) ! Saturation vapor pressure before PBL real(r8) :: ftem_prePBL(pcols,kte) ! Saturation vapor pressure before PBL real(r8) :: ftem_aftPBL(pcols,kte) ! Saturation vapor pressure after PBL real(r8) :: tem2( pcols,kte) ! Saturation specific humidity and RH real(r8) :: t_aftPBL( pcols,kte) ! Temperature after PBL diffusion real(r8) :: tten( pcols,kte) ! Temperature tendency by PBL diffusion real(r8) :: rhten( pcols,kte) ! RH tendency by PBL diffusion real(r8) :: qv_aft_PBL( pcols,kte) ! qv after PBL diffusion real(r8) :: ql_aft_PBL( pcols,kte) ! ql after PBL diffusion real(r8) :: qi_aft_PBL( pcols,kte) ! qi after PBL diffusion real(r8) :: s_aft_PBL( pcols,kte) ! s after PBL diffusion real(r8) :: u_aft_PBL( pcols,kte) ! u after PBL diffusion real(r8) :: v_aft_PBL( pcols,kte) ! v after PBL diffusion real(r8) :: qv_pro( pcols,kte) ! real(r8) :: ql_pro( pcols,kte) ! real(r8) :: qi_pro( pcols,kte) ! real(r8) :: s_pro( pcols,kte) ! real(r8) :: t_pro( pcols,kte) ! real(r8) :: u8( pcols,kte) ! x component of velocity in CAM's data structure (m/s) real(r8) :: v8( pcols,kte) ! y component of velocity in CAM's data structure (m/s) real(r8) :: t8( pcols,kte) ! real(r8) :: pmid8( pcols,kte) ! Pressure at the midpoints in CAM's data structure (Pa) real(r8) :: pmiddry8( pcols,kte) ! Dry Pressure at the midpoints in CAM's data structure (Pa) real(r8) :: zm8( pcols,kte) ! Height at the midpoints in CAM's data structure (m) real(r8) :: exner8( pcols,kte) ! exner function in CAM's data structure real(r8) :: s8( pcols,kte) ! Dry static energy (m2/s2) real(r8) :: rpdel8( pcols,kte) ! Inverse of pressure difference (1/Pa) real(r8) :: pdel8( pcols,kte) ! Pressure difference (Pa) real(r8) :: rpdeldry8( pcols,kte) ! Inverse of dry pressure difference (1/Pa) real(r8) :: pdeldry8( pcols,kte) ! dry pressure difference (1/Pa) REAL(r8) :: stnd( pcols,kte) ! Heating rate (dry static energy tendency, W/kg) real(r8) :: tke8( pcols,kte+1) ! Turbulent kinetic energy [ m2/s2 ] real(r8) :: turbtype( pcols,kte+1) ! Turbulent interface types [ no unit ] real(r8) :: smaw( pcols,kte+1) ! Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units] ! This is 1 when neutral condition (Ri=0), 4.964 for maximum unstable case, and 0 when Ri > Ricrit=0.19. real(r8) :: cgs( pcols,kte+1) ! Counter-gradient star [ cg/flux ] real(r8) :: cgh( pcols,kte+1) ! Counter-gradient term for heat [ J/kg/m ] real(r8) :: kvh( pcols,kte+1) ! Eddy diffusivity for heat [ m2/s ] real(r8) :: kvm( pcols,kte+1) ! Eddy diffusivity for momentum [ m2/s ] real(r8) :: kvq( pcols,kte+1) ! Eddy diffusivity for constituents [ m2/s ] real(r8) :: kvh_in( pcols,kte+1) ! kvh from previous timestep [ m2/s ] real(r8) :: kvm_in( pcols,kte+1) ! kvm from previous timestep [ m2/s ] real(r8) :: bprod( pcols,kte+1) ! Buoyancy production of tke [ m2/s3 ] real(r8) :: sprod( pcols,kte+1) ! Shear production of tke [ m2/s3 ] real(r8) :: sfi( pcols,kte+1) ! Saturation fraction at interfaces [ fraction ] real(r8) :: pint8( pcols,kte+1) ! Pressure at the interfaces in CAM's data structure (Pa) real(r8) :: pintdry8( pcols,kte+1) ! Dry pressure at the interfaces in CAM's data structure (Pa) real(r8) :: zi8( pcols,kte+1) ! Height at the interfacesin CAM's data structure (m) real(r8) :: cloud( pcols,kte,pcnst) ! Holder for cloud water and ice (q in cam) real(r8) :: cloudtnd( pcols,kte,pcnst) ! Holder for cloud tendencies (ptend_loc%q in cam) real(r8) :: wind_tends(pcols,kte,2) ! Wind component tendencies (m/s2) real(r8) :: cflx(pcols,pcnst) ! Surface constituent flux [ kg/m2/s ] #ifdef MODAL_AERO real(r8) :: tmp1(pcols) ! Temporary storage integer :: l, lspec #endif !! ----------------------- ! !! Main Computation Begins ! !! ----------------------- ! do icnst = 1 , pcnst !Setting curface fluxes for all constituents to be zero. Later, cflx(:,1) is populated by the water vapour surface flux cflx(:pcols,icnst) = 0.0_r8 enddo is_first_step = .false. if(itimestep == 1) then is_first_step = .true. endif ncol = pcols ztodt = DT rztodt = 1.0_r8 / ztodt !Few definitions in this subroutine require that ncol==1 if(ncol .NE. 1) then call wrf_error_fatal('Number of CAM Columns (NCOL) in CAMUWPBL scheme must be 1') endif !Initialize errstring errstring = '' !-------------------------------------------------------------------------------------! !Map CAM variables to the corresponding WRF variables ! !Loop over the points in the tile and treat them each as a CAM Chunk ! !-------------------------------------------------------------------------------------! itsp1 = its - 1 itile_len = ite - itsp1 do j = jts , jte do i = its , ite lchnk = (j - jts) * itile_len + (i - itsp1) !1-D index location from a 2-D tile phis(1) = ht(i,j) * gravit !Used for computing dry static energy !Flip vertically quantities computed at the mid points ktep1 = kte + 1 do k = kts,kte kflip = ktep1 - k !Loop is used as vector assignment may take more computational time do icnst = 1 , pcnst !Setting cloud and cloudtnd to values (obtained from module_cam_support) which should produce errors if used !1st,2nd,3rd and 5th constituents are used (see the assignments below) and diffused presently in this scheme cloud(1,kflip,icnst) = invalidVal cloudtnd(1,kflip,icnst) = invalidVal enddo u8( 1,kflip) = u_phy(i,k,j) ! X-component of velocity at the mid-points (m/s) [state%u in CAM] v8( 1,kflip) = v_phy(i,k,j) ! Y-component of velocity at the mid-points (m/s) [state%v in CAM] pmid8( 1,kflip) = p_phy(i,k,j) ! Pressure at the mid-points (Pa) [state%pmid in CAM] dp = p8w(i,k,j) - p8w(i,k+1,j) !Change in pressure (Pa) pdel8( 1,kflip) = dp rpdel8( 1,kflip) = 1.0_r8/dp ! Reciprocal of pressure difference (1/Pa) [state%rpdel in CAM] zm8( 1,kflip) = z(i,k,j)-ht(i,j) ! Height above the ground at the midpoints (m) [state%zm in CAM] t8( 1,kflip) = t_phy(i,k,j) ! Temprature at the mid points (K) [state%t in CAM] s8( 1,kflip) = cpair *t8(1,kflip) + gravit*zm8(1,kflip) + phis(1) ! Dry static energy (m2/s2) [state%s in CAM]-Formula tested in vertical_diffusion.F90 of CAM qrl8( 1,kflip) = rthratenlw(i,k,j)* cpair * dp ! Long Wave heating rate (W/kg*Pa)- Formula obtained from definition of qrlw(pcols,pver) in eddy_diff.F90 in CAM wsedl8( 1,kflip) = wsedl3d(i,k,j) ! Sedimentation velocity of stratiform liquid cloud droplet (m/s) !Following three formulas are obtained from ported CAM's ZM cumulus scheme !Values of 0 cause a crash in entropy multFrc = 1._r8/(1._r8 + qv_curr(i,k,j)) cloud( 1,kflip,1) = max( qv_curr(i,k,j) * multFrc, 1.0e-30_r8 ) !Specific humidity [state%q(:,:,1) in CAM] cloud( 1,kflip,ixcldliq) = qc_curr(i,k,j) * multFrc !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM] cloud( 1,kflip,ixcldice) = qi_curr(i,k,j) * multFrc !cloud ice [state%q(:,:,3) in CAM] cloud( 1,kflip,ixnumliq) = qnc_curr(i,k,j) * multFrc cloud( 1,kflip,ixnumice) = qni_curr(i,k,j) * multFrc exner8(1,kflip) = exner(i,k,j) !Exner function (no units) if(is_CAMMGMP_used) then cldn8( 1,kflip) = cldfra_old_mp(i,k,j) !Cloud fraction (no unit) else cldn8( 1,kflip) = cldfra(i,k,j) !Cloud fraction (no unit) endif !Following formula is obtained from physics_types.F90 of CAM (CESM101) pdeldry8(1,kflip) = pdel8(1,kflip)*(1._r8-cloud(1,kflip,1)) rpdeldry8(1,kflip) = 1._r8/pdeldry8(1,kflip) enddo do k = kts,kte+1 kflip = kte - k + 2 pint8( 1,kflip) = p8w( i,k,j) ! Pressure at interfaces [state%pint in CAM] zi8( 1,kflip) = z_at_w(i,k,j) -ht(i,j) ! Height at interfaces [state%zi in CAM] !Initialize Variables to zero, these are outputs from the "compute_eddy_diff" kvq(1,kflip) = 0.0_r8 ! Eddy diffusivity for constituents (m2/s) cgh(1,kflip) = 0.0_r8 ! Counter-gradient term for heat cgs(1,kflip) = 0.0_r8 ! Counter-gradient star (cg/flux) if( is_first_step ) then kvm3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for heat (m2/s) kvh3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for momentum (m2/s) endif kvh(1,kflip) = kvh3d(i,k,j) kvm(1,kflip) = kvm3d(i,k,j) end do !Compute pintdry8 and pmiddry8 !Following formula is obtained from physics_types.F90 of CAM (CESM101) pintdry8(1,1) = pint8(1,1) do k = 1, pver pintdry8(1,k+1) = pintdry8(1,k) + pdeldry8(1,k) pmiddry8(1,k) = (pintdry8(1,k+1)+ pintdry8(1,k))*0.5_r8 end do shflx( ncol) = hfx(i,j) ! Surface sensible heat flux (w/m2) !SGH and LANDFRAC are inputs for the compute_tms subroutine. Presently set to zero as do_tms is always false for now sgh( ncol) = 0.0_r8 ! Standard deviation of orography (m) landfrac(ncol) = 0.0_r8 ! Fraction (unitless) uMean = sqrt( u_phy(i,kts,j) * u_phy(i,kts,j) + v_phy(i,kts,j) * v_phy(i,kts,j) ) ! Mean velocity tauFac = rho(i,kts,j) * ustar(i,j) *ustar(i,j)/uMean taux(ncol) = -tauFac * u_phy(i,kts,j) ! x surface stress (N/m2) [Formulation obtained from CAM's BareGround.F90] tauy(ncol) = -tauFac * v_phy(i,kts,j) ! y surface stress (N/m2) !! Retrieve 'tauresx, tauresy' from from the last timestep if( is_first_step ) then tauresx2d(i,j) = 0._r8 tauresy2d(i,j) = 0._r8 endif tauresx(:ncol) = tauresx2d(i,j) tauresy(:ncol) = tauresy2d(i,j) !! All variables are modified by vertical diffusion !!------------------------------------------! !! Computation of turbulent mountain stress ! !!------------------------------------------! !! Consistent with the computation of 'normal' drag coefficient, we are using !! the raw input (u,v) to compute 'ksrftms', not the provisionally-marched 'u,v' !! within the iteration loop of the PBL scheme. if( do_tms ) then call compute_tms( pcols , pver , ncol , & u8 , v8 , t8 , pmid8 , & exner8 , zm8 , sgh , ksrftms , & tautmsx , tautmsy , landfrac ) !! Here, both 'taux, tautmsx' are explicit surface stresses. !! Note that this 'tautotx, tautoty' are different from the total stress !! that has been actually added into the atmosphere. This is because both !! taux and tautmsx are fully implicitly treated within compute_vdiff. !! However, 'tautotx, tautoty' are not used in the actual numerical !! computation in this module. tautotx(:ncol) = taux(:ncol) + tautmsx(:ncol) tautoty(:ncol) = tauy(:ncol) + tautmsy(:ncol) else ksrftms(:ncol) = 0.0_r8 tautotx(:ncol) = taux(:ncol) tautoty(:ncol) = tauy(:ncol) endif !-------------------------------------------------------------------------------------! !We are currenly using just water vapour flux at the surface, rest are set to zero !-------------------------------------------------------------------------------------! cflx(:pcols,1) = qfx(i,j) ! Surface constituent flux (kg/m2/s) !Following variables are initialized to zero, they are the output from the "compute_eddy_diff" call ustar8( :pcols) = 0.0_r8 ! Surface friction velocity (m/s) pblh( :pcols) = 0.0_r8 ! Planetary boundary layer height (m ) ipbl( :pcols) = 0.0_r8 ! If 1, PBL is CL, while if 0, PBL is STL. kpblh( :pcols) = 0.0_r8 ! Layer index containing PBL top within or at the base interface wstarPBL(:pcols) = 0.0_r8 ! Convective velocity within PBL (m/s) !!----------------------------------------------------------------------- ! !! Computation of eddy diffusivities - Select appropriate PBL scheme ! !!----------------------------------------------------------------------- ! select case (eddy_scheme) case ( 'diag_TKE' ) !! ---------------------------------------------------------------- ! !! At first time step, have eddy_diff.F90:caleddy() use kvh=kvm=kvf ! !! This has to be done in compute_eddy_diff after kvf is calculated ! !! ---------------------------------------------------------------- ! kvinit = .false. if(is_first_step) then kvinit = .true. endif !! ---------------------------------------------- ! !! Get LW radiative heating out of physics buffer ! !! ---------------------------------------------- ! !! Retrieve eddy diffusivities for heat and momentum from physics buffer !! from last timestep ( if first timestep, has been initialized by inidat.F90 ) kvm_in(:ncol,:) = kvm(:ncol,:) kvh_in(:ncol,:) = kvh(:ncol,:) call compute_eddy_diff( lchnk , pcols , pver , ncol , t8 , cloud(:,:,1) , ztodt , & cloud(:,:,2), cloud(:,:,3), s8 , rpdel8 , cldn8 , qrl8 , wsedl8 , zm8 , zi8 , & pmid8 , pint8 , u8 , v8 , taux , tauy , shflx , cflx(:,1), wstarent, & nturb , ustar8 , pblh , kvm_in , kvh_in , kvm , kvh , kvq , cgh , & cgs , tpert , qpert , wpert , tke8 , bprod , sprod , sfi , fqsatd , & kvinit , tauresx , tauresy, ksrftms, ipbl(:), kpblh(:), wstarPBL(:) , turbtype , smaw ) !! ----------------------------------------------- ! !! Store TKE in pbuf for use by shallow convection ! !! ----------------------------------------------- ! !! Store updated kvh, kvm in pbuf to use here on the next timestep do k = kts,kte+1 kflip = kte - k + 2 kvh3d(i,k,j) = kvh(1,kflip) kvm3d(i,k,j) = kvm(1,kflip) end do end select !!------------------------------------ ! !! Application of diffusivities ! !!------------------------------------ ! cloudtnd( :ncol,:,:) = cloud(:ncol,:,:) stnd( :ncol,: ) = s8( :ncol,: ) wind_tends(:ncol,:,1) = u8( :ncol,: ) wind_tends(:ncol,:,2) = v8( :ncol,: ) !!------------------------------------------------------ ! !! Write profile output before applying diffusion scheme ! !!------------------------------------------------------ ! sl_prePBL(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) & - ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice) qt_prePBL(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) & + cloudtnd(:ncol,:pver,ixcldice) call aqsat( t8, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver ) ftem_prePBL(:ncol,:) = cloud(:ncol,:,1)/ftem(:ncol,:)*100._r8 !! --------------------------------------------------------------------------------- ! !! Call the diffusivity solver and solve diffusion equation ! !! The final two arguments are optional function references to ! !! constituent-independent and constituent-dependent moleculuar diffusivity routines ! !! --------------------------------------------------------------------------------- ! !! Modification : We may need to output 'tautotx_im,tautoty_im' from below 'compute_vdiff' and !! separately print out as diagnostic output, because these are different from !! the explicit 'tautotx, tautoty' computed above. !! Note that the output 'tauresx,tauresy' from below subroutines are fully implicit ones. if( any(fieldlist_wet) ) then call compute_vdiff( lchnk, pcols, pver, pcnst, ncol, pmid8, pint8, rpdel8, t8, ztodt, & taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, qmincg, & fieldlist_wet, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, & tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, & compute_molec_diff, vd_lu_qdecomp) end if if( errstring .ne. '' ) then call wrf_error_fatal(errstring) endif if( any( fieldlist_dry ) ) then if( do_molec_diff ) then errstring = "Design flaw: dry vdiff not currently supported with molecular diffusion" call wrf_error_fatal(errstring) end if call compute_vdiff( lchnk, pcols, pver, pcnst, ncol, pmiddry8, pintdry8, rpdeldry8, t8, & ztodt, taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, & qmincg, fieldlist_dry, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, & tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, & compute_molec_diff, vd_lu_qdecomp) if( errstring .ne. '' ) call wrf_error_fatal(errstring) end if !! Store updated tauresx, tauresy to use here on the next timestep tauresx2d(i,j) = tauresx(ncol) tauresy2d(i,j) = tauresy(ncol) #ifdef MODAL_AERO !! Add the explicit surface fluxes to the lowest layer !! Modification : I should check whether this explicit adding is consistent with !! the treatment of other tracers. !The following code is commented out as the diffusion for the Aerosols and other species is handled by CAMMGMP and WRF_CHEM's dry_dep_driver subroutines !tmp1(:ncol) = ztodt * gravit * rpdel8(:ncol,pver) !do m = 1, ntot_amode ! l = numptr_amode(m) ! cloudtnd(:ncol,pver,l) = cloudtnd(:ncol,pver,l) + tmp1(:ncol) * cflx(:ncol,l) ! do lspec = 1, nspec_amode(m) ! l = lmassptr_amode(lspec,m) ! cloudtnd(:ncol,pver,l) = cloudtnd(:ncol,pver,l) + tmp1(:ncol) * cflx(:ncol,l) ! enddo !enddo #endif !! -------------------------------------------------------- ! !! Diagnostics and output writing after applying PBL scheme ! !! -------------------------------------------------------- ! sl(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) & - ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice) qt(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) & + cloudtnd(:ncol,:pver,ixcldice) !! --------------------------------------------------------------- ! !! Convert the new profiles into vertical diffusion tendencies. ! !! Convert KE dissipative heat change into "temperature" tendency. ! !! --------------------------------------------------------------- ! slten(:ncol,:) = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt qtten(:ncol,:) = ( qt(:ncol,:) - qt_prePBL(:ncol,:) ) * rztodt stnd(:ncol,:) = ( stnd(:ncol,:) - s8(:ncol,:) ) * rztodt wind_tends(:ncol,:,1) = ( wind_tends(:ncol,:,1) - u8(:ncol,:) ) * rztodt wind_tends(:ncol,:,2) = ( wind_tends(:ncol,:,2) - v8(:ncol,:) ) * rztodt cloudtnd(:ncol,:pver,:) = ( cloudtnd(:ncol,:pver,:) - cloud(:ncol,:pver,:) ) * rztodt !! ----------------------------------------------------------- ! !! In order to perform 'pseudo-conservative varible diffusion' ! !! perform the following two stages: ! !! ! !! I. Re-set (1) 'qvten' by 'qtten', and 'qlten = qiten = 0' ! !! (2) 'sten' by 'slten', and ! !! (3) 'qlten = qiten = 0' ! !! ! !! II. Apply 'positive_moisture' ! !! ! !! ----------------------------------------------------------- ! if( eddy_scheme .eq. 'diag_TKE' .and. do_pseudocon_diff ) then cloudtnd(:ncol,:pver,1) = qtten(:ncol,:pver) stnd(:ncol,:pver) = slten(:ncol,:pver) cloudtnd(:ncol,:pver,ixcldliq) = 0._r8 cloudtnd(:ncol,:pver,ixcldice) = 0._r8 cloudtnd(:ncol,:pver,ixnumliq) = 0._r8 cloudtnd(:ncol,:pver,ixnumice) = 0._r8 do ips = 1, ncol do k = 1, pver qv_pro(ips,k) = cloud(ips,k,1) + cloudtnd(ips,k,1) * ztodt ql_pro(ips,k) = cloud(ips,k,ixcldliq) + cloudtnd(ips,k,ixcldliq) * ztodt qi_pro(ips,k) = cloud(ips,k,ixcldice) + cloudtnd(ips,k,ixcldice) * ztodt s_pro(ips,k) = s8(ips,k) + stnd(ips,k) * ztodt t_pro(ips,k) = t8(ips,k) + (1._r8/cpair)*stnd(ips,k) * ztodt end do end do call positive_moisture( cpair, latvap, latvap+latice, ncol, pver, ztodt, qmin(1), qmin(2), qmin(3), & pdel8(:ncol,pver:1:-1), qv_pro(:ncol,pver:1:-1), ql_pro(:ncol,pver:1:-1), & qi_pro(:ncol,pver:1:-1), t_pro(:ncol,pver:1:-1), s_pro(:ncol,pver:1:-1), & cloudtnd(:ncol,pver:1:-1,1), cloudtnd(:ncol,pver:1:-1,ixcldliq), & cloudtnd(:ncol,pver:1:-1,ixcldice), stnd(:ncol,pver:1:-1) ) end if !! ----------------------------------------------------------------- ! !! Re-calculate diagnostic output variables after vertical diffusion ! !! ----------------------------------------------------------------- ! qv_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,1) + cloudtnd(:ncol,:pver,1) * ztodt ql_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldliq) + cloudtnd(:ncol,:pver,ixcldliq) * ztodt qi_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldice) + cloudtnd(:ncol,:pver,ixcldice) * ztodt s_aft_PBL(:ncol,:pver) = s8(:ncol,:pver) + stnd(:ncol,:pver) * ztodt t_aftPBL(:ncol,:pver) = ( s_aft_PBL(:ncol,:pver) - gravit*zm8(:ncol,:pver) ) / cpair u_aft_PBL(:ncol,:pver) = u8(:ncol,:pver) + wind_tends(:ncol,:pver,1) * ztodt v_aft_PBL(:ncol,:pver) = v8(:ncol,:pver) + wind_tends(:ncol,:pver,2) * ztodt call aqsat( t_aftPBL, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver ) ftem_aftPBL(:ncol,:pver) = qv_aft_PBL(:ncol,:pver) / ftem(:ncol,:pver) * 100._r8 tten(:ncol,:pver) = ( t_aftPBL(:ncol,:pver) - t8(:ncol,:pver) ) * rztodt rhten(:ncol,:pver) = ( ftem_aftPBL(:ncol,:pver) - ftem_prePBL(:ncol,:pver) ) * rztodt !Post processing of the output from CAM's parameterization do k=kts,kte kflip = kte-k+1 rublten(i,k,j) = wind_tends(1,kflip,1) rvblten(i,k,j) = wind_tends(1,kflip,2) rthblten(i,k,j) = stnd(1,kflip)/cpair/exner8(1,kflip) multFrc = 1._r8 + qv_curr(i,k,j) rqvblten(i,k,j) = cloudtnd(1,kflip,1 ) * multFrc * multFrc rqcblten(i,k,j) = cloudtnd(1,kflip,ixcldliq) * multFrc rqiblten(i,k,j) = cloudtnd(1,kflip,ixcldice) * multFrc !*Important* : ixnumliq is mixed in the dropmixnuc, therefore ixnumliq is NOT mixed here (ONLY if CAMMGMP is used for mp_physics) rqniblten(i,k,j) = cloudtnd(1,kflip,ixnumice) * multFrc !Diffusivity coefficients at the midpoints kp1 = k + 1 end do do k = kts,kte+1 kflip = kte - k + 2 tke_pbl(i,k,j) = tke8(1,kflip) !TKE at the interfaces turbtype3d(i,k,j) = turbtype(1,kflip) smaw3d(i,k,j) = smaw(1,kflip) end do kpbl2d(i,j) = kte - int(kpblh(1)) + 1 !int(kpblh(1)) pblh2d(i,j) = pblh( 1) tpert2d(i,j) = tpert(1) qpert2d(i,j) = qpert(1) wpert2d(i,j) = wpert(1) !End Post processing of the output from CAM enddo !loop of i enddo !loop of j return end subroutine camuwpbl !----------------------------------------------------------------------------------------- subroutine positive_moisture( cp, xlv, xls, ncol, mkx, dt, qvmin, qlmin, qimin, & dp, qv, ql, qi, t, s, qvten, qlten, qiten, sten ) !! ------------------------------------------------------------------------------- ! !! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, ! !! force them to be larger than minimum value by (1) condensating water vapor ! !! into liquid or ice, and (2) by transporting water vapor from the very lower ! !! layer. '2._r8' is multiplied to the minimum values for safety. ! !! Update final state variables and tendencies associated with this correction. ! !! If any condensation happens, update (s,t) too. ! !! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding ! !! input tendencies. ! !! Be careful the order of k : '1': near-surface layer, 'mkx' : top layer ! !! ------------------------------------------------------------------------------- ! !----------------------------------------------------------------------------------------- implicit none integer, intent(in) :: ncol, mkx real(r8), intent(in) :: cp, xlv, xls real(r8), intent(in) :: dt, qvmin, qlmin, qimin real(r8), intent(in) :: dp(ncol,mkx) real(r8), intent(inout) :: qv(ncol,mkx), ql(ncol,mkx), qi(ncol,mkx), t(ncol,mkx), s(ncol,mkx) real(r8), intent(inout) :: qvten(ncol,mkx), qlten(ncol,mkx), qiten(ncol,mkx), sten(ncol,mkx) integer i, k real(r8) dql, dqi, dqv, sum, aa, dum !! Modification : I should check whether this is exactly same as the one used in !! shallow convection and cloud macrophysics. do i = 1, ncol do k = mkx, 1, -1 ! From the top to the 1st (lowest) layer from the surface dql = max(0._r8,1._r8*qlmin-ql(i,k)) dqi = max(0._r8,1._r8*qimin-qi(i,k)) qlten(i,k) = qlten(i,k) + dql/dt qiten(i,k) = qiten(i,k) + dqi/dt qvten(i,k) = qvten(i,k) - (dql+dqi)/dt sten(i,k) = sten(i,k) + xlv * (dql/dt) + xls * (dqi/dt) ql(i,k) = ql(i,k) + dql qi(i,k) = qi(i,k) + dqi qv(i,k) = qv(i,k) - dql - dqi s(i,k) = s(i,k) + xlv * dql + xls * dqi t(i,k) = t(i,k) + (xlv * dql + xls * dqi)/cp dqv = max(0._r8,1._r8*qvmin-qv(i,k)) qvten(i,k) = qvten(i,k) + dqv/dt qv(i,k) = qv(i,k) + dqv if( k .ne. 1 ) then qv(i,k-1) = qv(i,k-1) - dqv*dp(i,k)/dp(i,k-1) qvten(i,k-1) = qvten(i,k-1) - dqv*dp(i,k)/dp(i,k-1)/dt endif qv(i,k) = max(qv(i,k),qvmin) ql(i,k) = max(ql(i,k),qlmin) qi(i,k) = max(qi(i,k),qimin) end do !! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally !! extracted from all the layers that has 'qv > 2*qvmin'. This fully !! preserves column moisture. if( dqv .gt. 1.e-20_r8 ) then sum = 0._r8 do k = 1, mkx if( qv(i,k) .gt. 2._r8*qvmin ) sum = sum + qv(i,k)*dp(i,k) enddo aa = dqv*dp(i,1)/max(1.e-20_r8,sum) if( aa .lt. 0.5_r8 ) then do k = 1, mkx if( qv(i,k) .gt. 2._r8*qvmin ) then dum = aa*qv(i,k) qv(i,k) = qv(i,k) - dum qvten(i,k) = qvten(i,k) - dum/dt endif enddo else write(iulog,*) 'Full positive_moisture is impossible in vertical_diffusion' call wrf_message(iulog) endif endif end do return end subroutine positive_moisture !----------------------------------------------------------------------------------------- subroutine camuwpblinit(rublten,rvblten,rthblten,rqvblten, & restart,tke_pbl,is_CAMMGMP_used, & ids,ide,jds,jde,kds,kde, & ims,ime,jms,jme,kms,kme, & its,ite,jts,jte,kts,kte) !!------------------------------------------------------------------! !! Initialization of time independent fields for vertical diffusion ! !! Calls initialization routines for subsidiary modules ! !!----------------------------------------------------------------- ! !This subroutine is based on vertical_diffusion_init of CAM. This subroutine !initializes variables for vertical diffusion subroutine calls. The layout !is kept similar to the original CAM subroutine to facillitate future adaptations. !Called by : physics_init.F !Calls: vd_register, cnst_get_ind, wrf_error_fatal,init_eddy_diff,init_tms,init_vdiff !----------------------------------------------------------------------------------------- use eddy_diff, only : init_eddy_diff use molec_diff, only : init_molec_diff use trb_mtn_stress, only : init_tms use diffusion_solver, only : init_vdiff, vdiff_select use constituents, only : cnst_get_ind, cnst_get_type_byind, cnst_name use module_cam_support, only : masterproc use module_model_constants, only : epsq2 #ifdef MODAL_AERO use modal_aero_data #endif implicit none !-------------------------------------------------------------------------------------! !Input and output variables ! !-------------------------------------------------------------------------------------! logical,intent(in) :: restart,is_CAMMGMP_used integer,intent(in) :: ids,ide,jds,jde,kds,kde integer,intent(in) :: ims,ime,jms,jme,kms,kme integer,intent(in) :: its,ite,jts,jte,kts,kte real,dimension(ims:ime,kms:kme,jms:jme),intent(inout) :: & rublten, rvblten, rthblten, rqvblten real,dimension(ims:ime,kms:kme,jms:jme),intent(out) :: TKE_PBL !-------------------------------------------------------------------------------------! !Local Variables ! !-------------------------------------------------------------------------------------! integer :: i,j,k,itf,jtf,ktf integer :: ntop_eddy !! Top interface level to which eddy vertical diffusion is applied ( = 1 ) integer :: nbot_eddy !! Bottom interface level to which eddy vertical diffusion is applied ( = pver ) integer :: ntop_molec !! Top interface level to which molecular vertical diffusion is applied ( = 1 ) integer :: nbot_molec !! Bottom interface level to which molecular vertical diffusion is applied character(128) :: errstring !! Error status for init_vdiff real(r8) :: hypm(kte) !! reference state midpoint pressures #ifdef MODAL_AERO integer :: m, l #endif jtf = min(jte,jde-1) ktf = min(kte,kde-1) itf = min(ite,ide-1) !Map CAM veritcal level variables pver = ktf - kts + 1 pverp = pver + 1 !Initialize flags and add constituents call vd_register() !! ----------------------------------------------------------------- ! !! Get indices of cloud liquid and ice within the constituents array ! !! ----------------------------------------------------------------- ! call cnst_get_ind( 'CLDLIQ', ixcldliq ) call cnst_get_ind( 'CLDICE', ixcldice ) if( microp_scheme .eq. 'MG' ) then call cnst_get_ind( 'NUMLIQ', ixnumliq ) call cnst_get_ind( 'NUMICE', ixnumice ) endif if (masterproc) then write(iulog,*)'Initializing vertical diffusion (vertical_diffusion_init)' call wrf_debug(1,iulog) end if !! ---------------------------------------------------------------------------------------- ! !! Initialize molecular diffusivity module ! !! Molecular diffusion turned on above ~60 km (50 Pa) if model top is above ~90 km (.1 Pa). ! !! Note that computing molecular diffusivities is a trivial expense, but constituent ! !! diffusivities depend on their molecular weights. Decomposing the diffusion matric ! !! for each constituent is a needless expense unless the diffusivity is significant. ! !! ---------------------------------------------------------------------------------------- ! ntop_molec = 1 !! Should always be 1 nbot_molec = 0 !! Should be set below about 70 km !! ---------------------------------- ! !! Initialize eddy diffusivity module ! !! ---------------------------------- ! ntop_eddy = 1 !! No reason not to make this 1, if > 1, must be <= nbot_molec nbot_eddy = pver !! Should always be pver if( masterproc ) write(iulog,fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY =', ntop_eddy, 'NBOT_EDDY =', nbot_eddy call wrf_debug(1,iulog) select case ( eddy_scheme ) case ( 'diag_TKE' ) if( masterproc ) write(iulog,*) 'vertical_diffusion_init: eddy_diffusivity scheme' call wrf_debug(1,iulog) if( masterproc ) write(iulog,*) 'UW Moist Turbulence Scheme by Bretherton and Park' call wrf_debug(1,iulog) !! Check compatibility of eddy and shallow scheme if( shallow_scheme .ne. 'UW' ) then write(iulog,*) 'ERROR: shallow convection scheme ', shallow_scheme,' is incompatible with eddy scheme ', eddy_scheme call wrf_message(iulog) call wrf_error_fatal( 'convect_shallow_init: shallow_scheme and eddy_scheme are incompatible' ) endif call init_eddy_diff( r8, pver, gravit, cpair, rair, zvir, latvap, latice, & ntop_eddy, nbot_eddy, hypm, karman ) if( masterproc ) write(iulog,*) 'vertical_diffusion: nturb, ntop_eddy, nbot_eddy ', nturb, ntop_eddy, nbot_eddy call wrf_debug(1,iulog) end select !!-------------------------------------------------------------------------------------! !! The vertical diffusion solver must operate ! !! over the full range of molecular and eddy diffusion ! !!-------------------------------------------------------------------------------------! ntop = min(ntop_molec,ntop_eddy) nbot = max(nbot_molec,nbot_eddy) !! ------------------------------------------- ! !! Initialize turbulent mountain stress module ! !! ------------------------------------------- ! if( do_tms ) then call init_tms( r8, tms_orocnst, tms_z0fac, karman, gravit, rair ) if (masterproc) then write(iulog,*)'Using turbulent mountain stress module' call wrf_message(iulog) write(iulog,*)' tms_orocnst = ',tms_orocnst call wrf_message(iulog) end if endif !! ---------------------------------- ! !! Initialize diffusion solver module ! !! ---------------------------------- ! call init_vdiff( r8, pcnst, rair, gravit, fieldlist_wet, fieldlist_dry, errstring ) if( errstring .ne. '' ) call wrf_error_fatal( errstring ) !!------------------------------------------------------------------------------------- !! Use fieldlist_wet to select the fields which will be diffused using moist mixing ratios ( all by default ) !! Use fieldlist_dry to select the fields which will be diffused using dry mixing ratios. !!------------------------------------------------------------------------------------- if( vdiff_select( fieldlist_wet, 'u' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'u' ) ) if( vdiff_select( fieldlist_wet, 'v' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'v' ) ) if( vdiff_select( fieldlist_wet, 's' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 's' ) ) #ifdef MODAL_AERO !Get index for the liquid number concentration (#/kg) call cnst_get_ind( 'NUMLIQ', ixndrop ) #endif do k = 1, pcnst #ifdef MODAL_AERO !! Do not diffuse droplet number - treated in dropmixnuc !BSINGH:01/31/2013: The following if condition is modified from its original form[which is:if( k == ixndrop ) cycle] !to allow the PBL scheme to mix liquid number when CAMMGMP scheme is not being used if(is_CAMMGMP_used) then!BSINGH:01/31/2013: Modified this if cond. see comment above if( k == ixndrop ) cycle endif !BSINGH-08/31/2012: In WRF, Physics init is called before chem init therefore !aerosol species are not yet added. Following condition, skip all other constituents !(which will be added in chem init) except first five (vapor, cld liq, cld ice, liq num and ice num) if( k > 5 ) cycle !The Following do-loop is commented out as Aerosols and other species are diffused in CAMMGMP and WRF_CHEM's dry_dep_driver subroutines !! Don't diffuse aerosol - treated in dropmixnuc !do m = 1, ntot_amode ! if( k == numptr_amode(m) ) cycle !! if( k == numptrcw_amode(m) ) go to 20 ! do l = 1, nspec_amode(m) ! if( k == lmassptr_amode(l,m) )cycle !! if( k == lmassptrcw_amode(l,m) ) go to 20 ! enddo !! if( k == lwaterptr_amode(m) ) go to 20 !enddo #endif !endif if( cnst_get_type_byind(k) .eq. 'wet' ) then if( vdiff_select( fieldlist_wet, 'q', k ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'q', k ) ) else !The DRY constituents are not diffused currently in PBL, they are diffused in WRF_CHEM's dry_dep_driver subroutine. Therefore the following line is commented out !if( vdiff_select( fieldlist_dry, 'q', k ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_dry, 'q', k ) ) endif enddo !Initialize the tendencies jtf=min0(jte,jde-1) ktf=min0(kte,kde-1) itf=min0(ite,ide-1) if(.not.restart)then do j = jts , jtf do k = kts , ktf do i = its , itf tke_pbl(i,k,j) = epsq2 rublten(i,k,j) = 0. rvblten(i,k,j) = 0. rthblten(i,k,j) = 0. rqvblten(i,k,j) = 0. enddo enddo enddo endif end subroutine camuwpblinit !----------------------------------------------------------------------------------------- subroutine vd_register() ! !This subroutine is based on the vd_register subroutine of CAM. !----------------------------------------------------------------------------------------- !Set flags for the PBL scheme (these flags are obtained using phys_getopts in CAM) microp_scheme = 'MG' !Used for adding constituents eddy_scheme = 'diag_TKE' !Used for calling eddy scheme !The following flag is deliberately set to UW for now. shallow_scheme = 'UW' !Eddy scheme is ONLY compaticle with 'UW' shallow scheme; do_tms = .false. !To include stresses due to orography tms_orocnst = 1._r8 !Orography constant end subroutine vd_register end module module_bl_camuwpbl_driver