!FDOC_TAG_GFDL module strat_cloud_legacy_mod !----------------------------------------------------------------------- ! NOTE: 10/7/10 -- RSH ! This is equivalent to the version of strat_cloud.F90 which existed ! before mods to include an option for double-moment microphysics (coded ! by m. salzmann) were added. ! ! This version is equivalent to that which was used in CM3 AR5 experiments. ! Treatment of internal diagnostics and the namelist within the module has ! been changed to be compatible with that used in the revised version of ! strat_cloud_mod, so that this code may be cleanly accessed from ! strat_cloud_mod while being in its own module, thus reducing compilation ! time. ! ! Addition of the double moment mods required some code changes which ! result in model answer changes when the code is optimized (-O2). However, ! when the atmos_phys library is compiled with -O0 and -fltconsistency, ! the results of an AM3 integration using the revised code (available in ! strat_cloud.F90) and this legacy version, may be seen to be identical. !----------------------------------------------------------------------- ! ! Stephen Klein ! !/"/>q ! ! Code to compute time tendencies of stratiform clouds and ! diagnoses ! rain and snow flux with prognostic scheme. ! ! ! ! ! ! The prognostic scheme returns the time tendencies of liquid, ! ice, and saturated volume fraction that are suspended in ! stratiform clouds. The scheme diagnoses the fluxes of rain ! and snow in saturated and unsaturated areas. ! ! The prognostic cloud scheme is responsible for determing ! cloud volume fractions, condensed water tendencies, and ! the stratiform precipitation rate. It includes processes ! for evaporation, condensation, deposition, and sublimation ! of cloud water, conversion of cloud water to precipitation, ! evaporation of falling precipitation, the bergeron-findeisan ! process, freezing of cloud liquid, accretion of cloud water ! by precipitation, and melting of falling precipitation. ! ! This scheme is based on the experience the author had ! at the ECMWF in 1997. The saturated volume fraction formalism ! and type of solution follow directly from the scheme of Tiedtke ! (1993): Monthly Weather Review, Volume 121, pages 3040-3061. ! The form of most of the microphysics follows Rotstayn , 1997: ! Quart. J. Roy. Met. Soc. vol 123, pages 1227-1282. The partial ! precipitation area formulism follows Jakob and Klein, 2000: ! Quart. J. Roy. Met. Soc. vol 126, pages 2525-2544. ! ! The statistical cloud scheme treatment, which is used as ! a replacement for the Tiedtke cloud fraction scheme, is based ! on a number of publications: Tompkins, A., 2002: J. Atmos. ! Sci., 59, 1917-1942, Klein et al., 2005: J. Geophys. Res., ! 110, D15S06, doi:10.1029/2004JD005017. ! ! ! ! native format of the restart file ! ! ! netcdf format of the restart file ! ! ! !The saturation volume fraction formalism comes from: ! !Tiedtke, M., 1993: Representation of clouds in large-scale models. ! Mon. Wea. Rev., 121, 3040-3061. ! ! ! !The form of most of the microphysics follows: ! !Rotstayn, L., 1997: A physically based scheme for the treatment of ! stratiform clouds and precipitation in large-scale models. I: ! Description and evaluation of microphysical processes. Quart. J. ! Roy. Met. Soc. 123, 1227-1282. ! ! ! ! ! ! ! !1. qmin should be chosen such that the range of {qmin, max(qa,ql !,qi)} is resolved by the precision of the numbers used. (default = ! 1.E-10) ! ! !2. Dmin will be MACHINE DEPENDENT and occur when ! ! !a. 1. -exp(-Dmin) = 0. instead of Dmin in the limit of very small ! Dmin ! !AND ! !b. 1. - exp(-D) < D for all D > Dmin ! ! ! use sat_vapor_pres_mod, only : compute_qs use beta_dist_mod, only: beta_dist_init, beta_dist_end, & incomplete_beta use aer_ccn_act_mod, only : aer_ccn_act_wpdf, aer_ccn_act_init use aer_in_act_mod, only : Jhete_dep use mpp_mod, only : mpp_clock_id, mpp_clock_begin, & mpp_clock_end, CLOCK_LOOP use fms_mod, only : file_exist, open_namelist_file,& & error_mesg, FATAL, NOTE, mpp_pe, mpp_root_pe, close_file,& & read_data, write_data, check_nml_error, write_version_number& &, stdlog, open_restart_file, open_ieee32_file, mpp_error use fms_io_mod, only : get_restart_io_mode, & register_restart_field, restart_file_type, & save_restart, restore_state, get_mosaic_tile_file use constants_mod, only : rdgas,rvgas,hlv,hlf,hls, cp_air& &,grav,tfreeze,dens_h2o use cloud_rad_mod, only : cloud_rad_init use diag_manager_mod, only : register_diag_field, send_data use time_manager_mod, only : time_type, get_date, get_time use cloud_generator_mod, only : do_cloud_generator, & cloud_generator_init, & compute_overlap_weighting use rad_utilities_mod, only : aerosol_type use strat_cloud_utilities_mod, & only: strat_cloud_utilities_init, & diag_id_type, diag_pt_type, & strat_nml_type, atmos_state_type, & particles_type implicit none private integer, private :: sc_loop, sc_pre_loop, sc_post_loop public strat_cloud_legacy, strat_cloud_legacy_init, & strat_cloud_legacy_end private aerosol_effects, ppm2m_sak, kmppm_sak, steepz_sak, & cloud_clear_xfer ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! GLOBAL STORAGE VARIABLES ! ! radturbten The sum of radiation and turbulent tendencies ! for each grid box. (K/sec) ! ! ! ------ data for cloud averaging code ------ ! real, allocatable, dimension (:,:,:) :: qlsum, qisum, cfsum integer, allocatable, dimension (:,:) :: nsum ! ! ------ constants used by the scheme ------- ! real, parameter :: d608 = (rvgas-rdgas) / rdgas real, parameter :: d622 = rdgas / rvgas real, parameter :: d378 = 1. - d622 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! DECLARE CONSTANTS AND SET DEFAULT VALUES FOR PARAMETERS OF ! THE SCHEME ! ! ! ! PHYSICAL CONSTANTS USED IN THE SCHEME ! ! ! constant definition unit ! ------------ ----------------------------- --------------- ! ! grav gravitational acceleration m/(s*s) ! ! hlv latent heat of vaporization J/kg condensate ! ! hlf latent heat of fusion J/kg condensate ! ! hls latent heat of sublimation J/kg condensate ! ! rdgas gas constant of dry air J/kg air/K ! ! rvgas gas constant of water vapor J/kg air/K ! ! cp_air specific heat of air at J/kg air/K ! constant pressure ! ! d622 rdgas divided by rvgas dimensionless ! ! d378 One minus d622 dimensionless ! ! tfreeze Triple point of water K ! ! dens_h2o Density of pure liquid kg/(m*m*m) ! ! ! ! ! PARAMETERS OF THE SCHEME ! ! ! parameter definition unit ! ------------ ----------------------------- --------------- ! ! ! rho_ice mass density of ice crystals kg/(m*m*m) ! ! ELI collection efficiency of dimensionless ! cloud liquid by falling ice ! ! ! do_average Average stratiform cloud properties ! before computing clouds used by radiation? ! ! overlap value of the overlap parameter ! from cloud rad ! overlap = 1 is maximum-random ! overlap = 2 is random ! real, parameter :: rho_ice = 100. real, parameter :: ELI = 0.7 logical :: do_average = .false. integer :: overlap = 2 ! !----------------------------------------------------------------------- !-------------------- diagnostics fields ------------------------------- !--- for netcdf restart type(restart_file_type), pointer, save :: Str_restart => NULL() type(restart_file_type), pointer, save :: Til_restart => NULL() logical :: in_different_file = .false. !---------------------------------------------------------- character(len=5) :: mod_name = 'strat' real :: missing_value = -999. real :: U00, rthresh, var_limit, sea_salt_scale, om_to_oc, N_land, & N_ocean, U_evap, eros_scale, eros_scale_c, eros_scale_t, & mc_thresh, diff_thresh, qmin, Dmin, efact, vfact, cfact, & iwc_crit, vfall_const2, vfall_exp2, qthalfwidth, N_min, & num_mass_ratio1, num_mass_ratio2 logical :: do_netcdf_restart, u00_profile, use_kk_auto, & use_online_aerosol, use_sub_seasalt, eros_choice, & super_choice, tracer_advec, do_old_snowmelt, retain_cm3_bug, do_pdf_clouds, & do_liq_num, do_dust_berg, pdf_org, do_ice_nucl_wpdf, debugo integer :: num_strat_pts, betaP, nsublevels, kmap, kord, & super_ice_opt, isamp, jsamp, ksamp character(len=64) :: microphys_scheme ! max_strat_pts maximum number of strat pts ! for instantaneous output ! integer, parameter:: max_strat_pts = 5 integer,dimension(2,max_strat_pts) :: strat_pts = 0 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! DECLARE VERSION NUMBER OF SCHEME ! Character(len=128) :: Version = '$Id: strat_cloud_legacy.F90,v 20.0 2013/12/13 23:22:13 fms Exp $' Character(len=128) :: Tagname = '$Name: tikal $' logical :: module_is_initialized = .false. integer, dimension(1) :: restart_versions = (/ 1 /) integer :: vers logical :: cloud_generator_on !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! The module contains the following subroutines: ! ! strat_cloud_legacy_init initializes module. ! ! strat_cloud_legacy calculations of the cloud scheme are performed ! here. ! ! strat_cloud_end closes out module. ! CONTAINS !####################################################################### subroutine strat_cloud_legacy_init (do_pdf_clouds) logical, intent(in) :: do_pdf_clouds if (module_is_initialized) return !----------------------------------------------------------------------- ! write version number to output file. !----------------------------------------------------------------------- call write_version_number(version, tagname) !----------------------------------------------------------------------- ! make sure needed modules have been initialized. !----------------------------------------------------------------------- call strat_cloud_utilities_init if (do_pdf_clouds) call beta_dist_init call cloud_generator_init cloud_generator_on = do_cloud_generator() !----------------------------------------------------------------------- ! mark this module initialized. !----------------------------------------------------------------------- module_is_initialized = .true. end subroutine strat_cloud_legacy_init !####################################################################### ! ! ! ! ! ! ! ! ! ! Time ! ! ! Indice of starting point in the longitude direction of the slab being passed to strat_cloud ! ! ! Indice of ending point in the longitude direction of the slab being passed ! ! ! Indice of starting point in the latitude direction of the slab being passed ! ! ! Indice of ending point in the latitude direction of the slab being passed ! ! ! Physics time step (sec) ! ! ! Pressure on model full levels (Pa) ! ! ! Pressure on model half levels (Pa) ! ! ! Sum of the tendencies of temperature from turbulence and radiation schemes (K/s) ! ! ! Temperature (K) ! ! ! Water vapor specific humidity (kg vapor/kg air) ! ! ! Grid-box mean liquid water specific humidity (kg liquid/kg air) ! ! ! Grid-box mean ice water specific humidity (kg ice/kg air) ! ! ! Cloud fraction (3d array and a prognostic variable) (fraction) ! ! ! Cloud droplet number (3d array and a prognostic variable) (#/kg air) ! ! ! Vertical pressure velocity (Pa/sec) ! ! ! Cumulus mass flux (defined positive as upward) (kg air/m2/sec) ! ! ! Vertical diffusion coefficient for temperature and tracer from vertical diffusion scheme (m2/sec) ! ! ! Fraction of surface that contains land (fraction) ! ! ! Change in temperature due to strat_cloud (K) ! ! ! Change in water vapor due to strat_cloud (kg vapor/kg air) ! ! ! Change in cloud liquid due to strat_cloud (kg liquid/kg air) ! ! ! Change in cloud ice due to strat_cloud (kg ice/kg air) ! ! ! Change in cloud fraction due to strat_cloud (fraction) ! ! ! Change in cloud droplet number due to strat_cloud (fraction) ! ! ! Surface rain fall over time step dtcloud (kg liquid/m2) ! ! ! Surface snow fall over time step dtcloud (kg ice/m2) ! ! ! 3D rain fall over time step dtcloud (kg liquid/m2) ! ! ! 3D snow fall over time step dtcloud (kg ice/m2) ! ! ! Ratio of large-scale specific humidity to specific humidity in ! environment outside convective system (from donner_deep) ! ! Will be equal to 1 for all normal AM2 operations (i.e. donner_deep is not activated) ! ! Note that index 1 is nearest ground ! ! ! ! The fraction of the grid box containing either cumulus cells ! or the mesoscale circulation (from donner_deep). ! ! Will be equal to 0 for all normal AM2 operations ! (i.e. donner_deep is not activated) ! ! Note that index 1 is nearest ground ! ! ! ! Optional input real array indicating the point is above the surface ! if equal to 1.0 and indicating the point is below the surface if ! equal to 0. ! ! Used only in eta vertical coordinate model. ! ! ! subroutine strat_cloud_legacy( Nml, & diag_id, diag_pt, n_diag_4d, n_diag_4d_kp1, & diag_4d, diag_4d_kp1, diag_3d, & Time,is,ie,js,je,dtcloud,pfull,phalf,radturbten2,& T,qv,ql,qi,qa,omega,Mc,diff_t,LAND, & ST,SQ,SL,SI,SA,f_snow_berg,rain3d,snow3d,snowclr3d,surfrain, & surfsnow,qrat,ahuco,limit_conv_cloud_frac,MASK, & qn, Aerosol, SN) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! ! VARIABLES ! ! ! ! ------ ! INPUT: ! ------ ! ! variable definition unit ! ------------ ----------------------------- --------------- ! ! Time time type variable ! ! is,ie starting and ending i indices ! for data window ! ! js,je starting and ending j indices ! for data window ! ! dtcloud time between this call and s ! the next call to strat_cloud ! ! pfull pressure at full model levels Pa ! IMPORTANT NOTE: p(j) 0) then strat_pts(:,1:) = Nml%strat_pts(:,1:) endif efact = Nml%efact vfact = Nml%vfact cfact = Nml%cfact do_old_snowmelt = Nml%do_old_snowmelt retain_cm3_bug = Nml%retain_cm3_bug do_pdf_clouds = Nml%do_pdf_clouds betaP = Nml%betaP iwc_crit = Nml%iwc_crit vfall_const2 = Nml%vfall_const2 vfall_exp2 = Nml%vfall_exp2 qthalfwidth = Nml%qthalfwidth nsublevels = Nml%nsublevels kmap = Nml%kmap kord = Nml%kord do_liq_num = Nml%do_liq_num do_dust_berg = Nml%do_dust_berg N_min = Nml%N_min num_mass_ratio1 = Nml%num_mass_ratio1 num_mass_ratio2 = Nml%num_mass_ratio2 microphys_scheme = Nml%microphys_scheme super_ice_opt = Nml%super_ice_opt pdf_org = Nml%pdf_org do_ice_nucl_wpdf = Nml%do_ice_nucl_wpdf debugo = Nml%debugo isamp = Nml%isamp jsamp = Nml%jsamp ksamp = Nml%ksamp !----------------------------------------------------------------------- ! ! if (.not.module_is_initialized) then call error_mesg('strat_cloud_legacy','strat_cloud_legacy is not initialized',FATAL) endif if (diag_id%ql_wt > 0) then diag_4d(:,:,:,diag_pt%ql_wt) = ql endif !----------------------------------------------------------------------- ! ! initialize select variables to zero. The variables reset ! are: ! ! (1) changes of prognostic variables ! ! (2) variables dealing with the rain/snow fluxes. ! ! (3) qa_mean of the level above the top level. ! ! (4) diagnostic output fields ST = 0. SQ = 0. SL = 0. SI = 0. SA = 0. if (present(SN)) SN = 0. rain_cld = 0. rain_clr = 0. a_rain_cld = 0. a_rain_clr = 0. snow_cld = 0. snow_clr = 0. a_snow_clr = 0. a_snow_cld = 0. qa_mean_lst= 0. dcond_ls = 0. dcond_ls_ice = 0. qcg = 0. qcg_ice = 0. !----------------------------------------------------------------------- ! ! Determine dimensions of slab idim = SIZE(T,1) jdim = SIZE(T,2) kdim = SIZE(T,3) !----------------------------------------------------------------------- ! ! compute inverse time step inv_dtcloud = 1.0 / dtcloud !----------------------------------------------------------------------- ! ! Calculate saturation specific humidity and its temperature ! derivative, and thermal conductivity plus vapor diffusivity ! factor. ! ! These are calculated according to the formulas: ! ! (1) qs = d622*esat/ [pfull - (1.-d622)*esat] ! ! (2) dqsdT = d622*pfull*(desat/dT)/[pfull-(1.-d622)*esat]**2. ! ! (3) gamma = (L/cp) * dqsdT ! ! where d622 = rdgas/rvgas; esat = saturation vapor pressure; ! and desat/dT is the temperature derivative of esat. ! Note that in the calculation of gamma, ! ! { hlv for T > tfreeze } ! L = { 0.05*(T-tfreeze+20.)*hlv + 0.05*(tfreeze-T)*hls } ! { for tfreeze-20.< T < tfreeze} ! { hls for T < tfreeze-20. } ! ! This linear form is chosen because at tfreeze-20. es = esi, and ! at tfreeze, es = esl, with linear interpolation in between. ! ! The conductivity/diffusivity factor, A_plus_B is given by: ! ! (4) A_plus_B = { (hlv/Ka/T)*((hlv/rvgas/T)-1.) } + ! ! { (rvgas*T/chi*esat) } ! ! where Ka is the thermal conductivity of air = 0.024 J/m/s/K ! and chi is the diffusitivy of water vapor in air which is ! given by ! ! (5) chi = 2.21 E-05 (m*m)/s * (1.E+05)/pfull ! ! where p is the pressure in Pascals. ! ! ! Note that qs, dqsdT, and gamma do not have their proper values ! until all of the following code has been executed. That ! is qs and dqsdT are used to store intermediary results ! in forming the full solution. !calculate water saturated vapor pressure from table !and store temporarily in the variable gamma !calculate qs and dqsdT if (do_pdf_clouds) then call compute_qs( T-((hlv*ql+hls*qi)/cp_air), pfull, qs, & dqsdT=dqsdT, esat=gamma ) else call compute_qs( T, pfull, qs, dqsdT=dqsdT, esat=gamma ) end if !compute A_plus_B A_plus_B = ( (hlv/0.024/T) * ((hlv/rvgas/T)-1.) ) + & (rvgas*T*pfull/2.21/gamma) !calculate gamma if (do_pdf_clouds) then gamma = dqsdT *(min(1.,max(0.,0.05*(T-((hlv*ql+hls*qi)/cp_air) & -tfreeze+20.)))*hlv + & min(1.,max(0.,0.05*(tfreeze -T+((hlv*ql+hls*qi)& /cp_air) )))*hls)/cp_air else gamma = dqsdT *(min(1.,max(0.,0.05*(T-tfreeze+20.)))*hlv + & min(1.,max(0.,0.05*(tfreeze -T )))*hls)/cp_air end if !----------------------------------------------------------------------- ! ! Calculate air density ! airdens = pfull / (rdgas * T * (1. + d608*qv - ql - qi) ) airdens = pfull / (rdgas * T * (1. - ql - qi) ) where (qrat .gt. 0.) airdens = pfull / (rdgas * T *(1.+(d608*qv/qrat)-ql-qi) ) end where !----------------------------------------------------------------------- ! ! Assign cloud droplet number based on land or ocean point. N = N_land*LAND + N_ocean*(1.-LAND) !--------------------------------------------------------------------- ! call aerosol_effects to include impact of aerosols on the cloud ! droplet number and the bergeron process, if these effects activated. !--------------------------------------------------------------------- if (do_liq_num .or. do_dust_berg) then call aerosol_effects (is, js, Time, phalf, airdens, T, & diag_4d, diag_id, diag_pt, & concen_dust_sub, totalmass1, Aerosol, mask) endif !----------------------------------------------------------------------- ! ! Is a sub-vertical grid scale distribution going to be neededed? ! If yes, then do ppm fits ! if (do_pdf_clouds) then !initialize quantities do j = 1, kdim delp(:,:,j) = phalf(:,:,j+1)-phalf(:,:,j) enddo qta4(1,:,:,:) = max(qmin,qv+ql+qi) qtqsa4(1,:,:,:) = qta4(1,:,:,:)-qs if (nsublevels.gt.1) then do id = 1,idim call ppm2m_sak(qta4(:,id,:,:),delp(id,:,:),kdim,kmap,1,jdim,0,kord) call ppm2m_sak(qtqsa4(:,id,:,:),delp(id,:,:),kdim,kmap,1,jdim,0,kord) enddo else qta4(2,:,:,:) = qta4(1,:,:,:) qta4(3,:,:,:) = qta4(1,:,:,:) qta4(4,:,:,:) = 0. qtqsa4(2,:,:,:) = qtqsa4(1,:,:,:) qtqsa4(3,:,:,:) = qtqsa4(1,:,:,:) qtqsa4(4,:,:,:) = 0. end if end if !end of do_pdf_clouds section !rab - statements moved outside of large loop ! they have no bearing on the overall cloud physics ! they are merely diagnostic values ! if (.not. do_pdf_clouds) then if (max(diag_id%qadt_fill,diag_id%qa_fill_col) > 0) then where (qa .le. qmin) ! qadt_fill = -qa * inv_dtcloud diag_4d(:,:,:,diag_pt%qadt_fill) = -qa * inv_dtcloud endwhere endif if (max(diag_id%qidt_fill,diag_id%qi_fill_col) > 0) then where (qi.le.qmin .or. qa.le.qmin) ! qidt_fill = -qi * inv_dtcloud diag_4d(:,:,:,diag_pt%qidt_fill) = -qi * inv_dtcloud endwhere endif if (.not. do_liq_num) then !a N3D = 0. diag_4d(:,:,:,diag_pt%droplets) = 0. if (max(diag_id%qldt_fill,diag_id%ql_fill_col) > 0) then where (ql .le. qmin .or. qa .le. qmin) ! qldt_fill = -ql * inv_dtcloud diag_4d(:,:,:,diag_pt%qldt_fill) = -ql * inv_dtcloud endwhere endif else !a N3D = qn*airdens*1.e-6 diag_4d(:,:,:,diag_pt%droplets) = qn(:,:,:)*airdens(:,:,:)*1.e-6 ! if (diag_id%debug1_3d > 0) debug1 = min(qa,1.) if (diag_id%cf_liq_init > 0) diag_4d(:,:,:,diag_pt%cf_liq_init ) = min(qa,1.) do j=1,kdim do k=1,jdim do i=1,idim if (ql(i,k,j).le.qmin .or. qa(i,k,j).le.qmin .or. qn(i,k,j).le.qmin) then !a N3D(i,k,j) = 0. diag_4d(i,k,j,diag_pt%droplets) = 0. ! if (max(diag_id%qldt_fill,diag_id%ql_fill_col) > 0) qldt_fill(i,k,j) = -ql(i,k,j) * inv_dtcloud if (max(diag_id%qldt_fill,diag_id%ql_fill_col) > 0) diag_4d(i,k,j,diag_pt%qldt_fill) = -ql(i,k,j) * inv_dtcloud ! if (max(diag_id%qndt_fill,diag_id%qn_fill_col) > 0) qndt_fill(i,k,j) = -qn(i,k,j) * inv_dtcloud if (max(diag_id%qndt_fill,diag_id%qn_fill_col, & diag_id%qldt_fill,diag_id%ql_fill_col) > 0) & diag_4d(i,k,j,diag_pt%qndt_fill) = -qn(i,k,j) * inv_dtcloud ! if (diag_id%debug1_3d > 0) debug1(i,k,j) = 0. if (diag_id%cf_liq_init > 0) diag_4d(i,k,j,diag_pt%cf_liq_init ) = 0. endif enddo enddo enddo endif else if (max(diag_id%qidt_fill,diag_id%qi_fill_col) > 0) then where (qi .le. qmin) ! qidt_fill = -qi * inv_dtcloud diag_4d(:,:,:,diag_pt%qidt_fill) = -qi * inv_dtcloud endwhere endif if (.not. do_liq_num) then !a N3D = 0. diag_4d(:,:,:,diag_pt%droplets) = 0. if (max(diag_id%qldt_fill,diag_id%ql_fill_col) > 0) then where (ql .le. qmin) ! qldt_fill = -ql * inv_dtcloud diag_4d(:,:,:,diag_pt%qldt_fill) = -ql * inv_dtcloud endwhere endif else !a N3D = qn*airdens*1.e-6 diag_4d(:,:,:,diag_pt%droplets) = qn(:,:,:)*airdens(:,:,:)*1.e-6 ! if (diag_id%debug1_3d > 0) debug1 = min(qa,1.) if (diag_id%cf_liq_init > 0) diag_4d(:,:,:,diag_pt%cf_liq_init ) = min(qa,1.) do j=1,kdim do k=1,jdim do i=1,idim if (ql(i,k,j).le.qmin .or. qn(i,k,j).le.qmin) then !a N3D(i,k,j) = 0. diag_4d(i,k,j,diag_pt%droplets) = 0. ! if (max(diag_id%qldt_fill,diag_id%ql_fill_col) > 0) qldt_fill(i,k,j) = -ql(i,k,j) * inv_dtcloud if (max(diag_id%qldt_fill,diag_id%ql_fill_col) > 0) diag_4d(i,k,j,diag_pt%qldt_fill) = -ql(i,k,j) * inv_dtcloud ! Should this just be diag_id%qndt_fill + diag_id%qn_fill_col ? !ORIGINAL: if (diag_id%qldt_fill > 0) & if (max(diag_id%qndt_fill,diag_id%qn_fill_col, & diag_id%qldt_fill,diag_id%ql_fill_col) > 0) & ! qndt_fill(i,k,j) = -qn(i,k,j) * inv_dtcloud diag_4d(i,k,j,diag_pt%qndt_fill) = -qn(i,k,j) * inv_dtcloud ! if (diag_id%debug1_3d > 0) debug1(i,k,j) = 0. if (diag_id%cf_liq_init > 0) diag_4d(i,k,j,diag_pt%cf_liq_init ) = 0. endif enddo enddo enddo endif endif !rab - end of statements moved outside large loop !a diag_4d(:,:,:, diag_pt%droplets) = N3D ! diag_4d(:,:,:, diag_pt%droplets_wtd) = N3D*ql call mpp_clock_end(sc_pre_loop) call mpp_clock_begin(sc_loop) !----------------------------------------------------------------------- ! ! Enter the large loop over vertical levels. Level 1 is the top ! level of the column and level kdim is the bottom of the model. ! If MASK is present, each column may not have kdim valid levels. ! rain3d = 0. snow3d = 0. snowclr3d = 0. DO j = 1, kdim sum_freeze = 0. sum_rime = 0. sum_berg = 0. f_snow_berg(:,:,j) = 0. !----------------------------------------------------------------------- ! ! Calculate pressure thickness of level and relative humidity !calculate difference in pressure across the grid box divided !by gravity deltpg = (phalf(:,:,j+1)-phalf(:,:,j))/grav !calculate GRID box mean relative humidity ! U = min(max(0.,qv(:,:,j)/qs(:,:,j)),1.) U = 0. where (qrat(:,:,j) .gt. 0.) U = min(max(0.,(qv(:,:,j)/(qrat(:,:,j)*qs(:,:,j)))),1.) end where !----------------------------------------------------------------------- ! ! Account for the fact that other processes may have created ! negative tracer or extremely small values of tracer fields. ! The general reason for the extremely small values of the ! tracer fields is due to vertical diffusion, advection of ! condensate or cumulus induced subsidence (also a form of ! advection) of condensate. ! ! In this step any values of the prognostic variables which are ! less than qmin are reset to zero, while conserving total ! moisture. ! ! Note that this is done slightly different for the Tiedtke ! cloud fraction than it is for pdf clouds. In the former, ! the filling requires that cloud liquid, cloud ice, and ! cloud fraction are greater than qmin. For PDF clouds, ! cloud fraction need not be considered since it is diagnosed ! below from the PDF clouds ! if (.not. do_pdf_clouds) then do k=1,jdim do i=1,idim qa_upd(i,k) = qa(i,k,j) qi_upd(i,k) = qi(i,k,j) ql_upd(i,k) = ql(i,k,j) if (qa(i,k,j) .le. qmin) then SA(i,k,j) = SA(i,k,j) - qa(i,k,j) qa_upd(i,k) = 0. end if ! Correct for qa > RH, which is not permitted under the ! assumption that the cloudy air is saturated and the temperature ! inside and outside of the cloud are about the same. if (qa_upd(i,k) .gt. U(i,k)) then ! if (max(diag_id%qadt_rhred,diag_id%qa_rhred_col) > 0 ) qadt_rhred(i,k,j) = (qa_upd(i,k)-U(i,k)) * inv_dtcloud if (max(diag_id%qadt_rhred,diag_id%qa_rhred_col) > 0 ) diag_4d(i,k,j,diag_pt%qadt_rhred) = (qa_upd(i,k)-U(i,k)) * inv_dtcloud SA(i,k,j) = SA(i,k,j) + U(i,k) - qa_upd(i,k) qa_upd(i,k) = U(i,k) end if if (.not. do_liq_num) then if (ql(i,k,j) .le. qmin .or. qa(i,k,j) .le. qmin) then SL(i,k,j) = SL(i,k,j) - ql(i,k,j) SQ(i,k,j) = SQ(i,k,j) + ql(i,k,j) ST(i,k,j) = ST(i,k,j) - hlv*ql(i,k,j)/cp_air ql_upd(i,k) = 0. end if else qn_upd(i,k) = qn(i,k,j) if (ql(i,k,j) .le. qmin .or. qa(i,k,j) .le. qmin .or. qn(i,k,j) .le. qmin) then SL(i,k,j) = SL(i,k,j) - ql(i,k,j) SQ(i,k,j) = SQ(i,k,j) + ql(i,k,j) ST(i,k,j) = ST(i,k,j) - hlv*ql(i,k,j)/cp_air SN(i,k,j) = SN(i,k,j) - qn(i,k,j) ql_upd(i,k) = 0. qn_upd(i,k) = 0. endif endif if (qi(i,k,j) .le. qmin .or. qa(i,k,j) .le. qmin) then SI(i,k,j) = SI(i,k,j) - qi(i,k,j) SQ(i,k,j) = SQ(i,k,j) + qi(i,k,j) ST(i,k,j) = ST(i,k,j) - hls*qi(i,k,j)/cp_air qi_upd(i,k) = 0. endif enddo enddo else do k=1,jdim do i=1,idim ql_upd(i,k) = ql(i,k,j) qi_upd(i,k) = qi(i,k,j) if (.not. do_liq_num) then if (ql(i,k,j) .le. qmin) then SL(i,k,j) = SL(i,k,j) - ql(i,k,j) SQ(i,k,j) = SQ(i,k,j) + ql(i,k,j) ST(i,k,j) = ST(i,k,j) - hlv*ql(i,k,j)/cp_air ql_upd(i,k) = 0. endif else qn_upd(i,k) = qn(i,k,j) if (ql(i,k,j) .le. qmin .or. qn(i,k,j) .le. qmin) then SL(i,k,j) = SL(i,k,j) - ql(i,k,j) SQ(i,k,j) = SQ(i,k,j) + ql(i,k,j) ST(i,k,j) = ST(i,k,j) - hlv*ql(i,k,j)/cp_air SN(i,k,j) = SN(i,k,j) - qn(i,k,j) ql_upd(i,k) = 0. qn_upd(i,k) = 0. endif endif if (qi(i,k,j) .le. qmin) then SI(i,k,j) = SI(i,k,j) - qi(i,k,j) SQ(i,k,j) = SQ(i,k,j) + qi(i,k,j) ST(i,k,j) = ST(i,k,j) - hls*qi(i,k,j)/cp_air qi_upd(i,k) = 0. endif end do end do end if !for do_pdf_clouds !----------------------------------------------------------------------! ! ! ! ! ! ! ! NON-CONVECTIVE CONDENSATION ! ! ! ! ! ! ! ! METHOD 1 ! ! ! ! ! ! TIEDTKE (1993) CLOUD FRACTION ! ! ! ! ANALYTIC INTEGRATION OF SATURATED VOLUME FRACTION EQUATION ! ! ! ! ! ! ! Do non-convective condensation following Tiedtke, pages 3044-5. ! In this formulation stratiform clouds are only formed/destroyed ! when there is upward or downward motion to support/destroy it. ! ! The first step is to compute the change in qs due to large- ! scale processes, dqs_ls. In Tiedtke, it has contributions from ! large-scale uplift, convection induced compensating subsidence, ! turbulence cooling and radiative cooling. dqs_ls has the form: ! ! (((omega+ grav*Mc)/airdens/cp)+radturbten)*dqsdT*dtcloud ! (6) dqs_ls= -------------------------------------------------------- ! 1. + ( qa + (da_ls/2.) ) * gamma ! ! Here da_ls is the increase in cloud fraction due to non- ! convective processes. Because this increase is also a function ! of dqs_ls, a quadratic equation must be solved for dqs_ls in ! the case that da_ls is not equal to zero. ! ! Note that if the PDF cloud scheme is active the Tiedtke large- ! scale condensation is bypassed. if (.not.do_pdf_clouds) then do k=1,jdim do i=1,idim dqs_ls(i,k) =(((omega(i,k,j)+grav*Mc(i,k,j))/airdens(i,k,j)/cp_air)+& radturbten2(i,k,j))*dtcloud*dqsdT(i,k,j) !compute pressure dependent U00 following ECMWF formula if !desired U00p(i,k) = U00 if (u00_profile) then if (pfull(i,k,j) .gt. 0.8*phalf(i,k,KDIM+1)) then U00p(i,k) = U00 + (1.-U00)* & (((pfull(i,k,j)-(0.8*phalf(i,k,KDIM+1))) & / (0.2*phalf(i,k,KDIM+1)) )**2.) endif end if !ljd ! modify u00p to account for humidity in convective system ! See "Tiedtke u00 adjustment" notes, 10/22/02 !ljd u00p(i,k)=u00p(i,k)+(1.-u00p(i,k))*ahuco(i,k,j) if (dqs_ls(i,k).le.0. .and. U(i,k).ge.U00p(i,k) .and. qa_upd(i,k).lt.1.) then tmp1s = sqrt( (1.+qa_upd(i,k)*gamma(i,k,j))**2. - (1.-qa_upd(i,k)) * & (1.-qa_upd(i,k))*gamma(i,k,j)*dqs_ls(i,k)/qs(i,k,j)/ & max(1.-U(i,k),qmin) ) - (1.+qa_upd(i,k)*gamma(i,k,j)) tmp1s = -1. * tmp1s / ((1.-qa_upd(i,k))*(1.-qa_upd(i,k))*gamma(i,k,j)/& qs(i,k,j)/max(1.-U(i,k),qmin)/2.) dqs_ls(i,k) = min(tmp1s,dqs_ls(i,k)/(1.+0.5*(1.+qa_upd(i,k))*gamma(i,k,j))) else dqs_ls(i,k) = dqs_ls(i,k)/(1.+qa_upd(i,k)*gamma(i,k,j)) endif ! The next step is to compute the change in saturated volume ! fraction due to non-convective condensation, da_ls. This ! occurs in two conditions: ! ! (a) dqs_ls < 0. and U00 < U < 1., where U00 is the threshold ! relative humidity for non-convective condensation. Note ! that if U is greater than or equal to 1., ideally qa = 1, ! and da_ls = 0. However this may not be the case for ! numerical reasons so this must be assured after analytic ! integration of the qa equation. ! ! For these cases the change in saturated volume fraction is: ! ! (7) da_ls = - (1.-qa)*(1.-qa)*dqs_ls/2./qs/(1.-U) ! ! This formula arises from the assumption that vapor is uni- ! formly distributed in the range [qv_clr - (qs - qv_clr),qs] ! where qv_clr is the amount of vapor in the unsaturated ! volume and is given from the following equation: ! ! (8) qv = qa * qs + (1.-qa) * qv_clr ! ! Implicit in equation (7) is the following assumption: ! As qsat changes, the distribution of qv+ql+qi ! remains constant. That is as qsat rises, portions where ! qv+ql+qi > qsat+dqsat remain saturated. This can only ! occur if it is assumed that ql+qi evaporate-sublimate or ! condense-deposit to keep qv = qsat. ! ! (b) dqs_ls > 0. Ideally some portion of the cloud should ! evaporate however this is not accounted for at present. ! !compute formula for da_ls ! where (dqs_ls .le. 0. .and. U .ge. U00p) ! da_ls = -0.5 * (1.-qa_upd) * (1.-qa_upd) * dqs_ls / & ! qs(:,:,j) / max(1.-U,qmin) ! elsewhere ! da_ls = 0. ! end where da_ls(i,k) = 0. if ((dqs_ls(i,k).le.0. .and. U(i,k).ge.U00p(i,k)) .and. & (qa_upd(i,k)+ahuco(i,k,j).le.1.)) then da_ls(i,k) = -0.5 * (1.-qa_upd(i,k)-ahuco(i,k,j)) * (1.-qa_upd(i,k)- & ahuco(i,k,j)) * dqs_ls(i,k)/ qs(i,k,j) / max(1.-U(i,k),qmin) endif enddo enddo ! Turbulent erosion of clouds ! ! As in Tiedtke (1993) this is calculated using the eros_scale ! parameter as: ! ! (9) dql/dt = - qa * eros_scale * (qs - qv) * (ql/ ql+qi ) ! ! (10) dqi/dt = - qa * eros_scale * (qs - qv) * (qi/ ql+qi ) ! ! (11) dqa/dt = - qa * eros_scale * (qs - qv) * (qa/ ql+qi ) ! ! for which the erosion sink term (B in equation 13) is ! ! (12) B = qa * eros_scale * (qs - qv) / (ql + qi) ! ! ! Theory for eros_scale ! ! If eros_choice equals false, then a single erosion time scale ! is used in all conditions (eros_scale). If eros_choice equals ! true then it is assumed that the timescale for turbulent ! evaporation is a function of the conditions in the grid box. ! Specifically, if the flow is highly turbulent then the scale is ! short, and eros_scale is large. Likewise if convection is ! occurring, then it is assumed that the erosion term is larger ! than backround conditions. ! ! Here are the typical values for the timescales and the ! switches used: ! ! Mixing type eros_scale (sec-1) Indicator ! ---------------- ------------------ -------------------- ! ! Background 1.e-06 always present ! Convective layers 5.e-06 Mc > Mc_thresh ! Turbulent layers 5.e-05 diff_t > diff_thresh ! do k=1,jdim do i=1,idim !Background erosion scale tmp2s = eros_scale !Do enhanced erosion in convective or turbulent layers? ! ! IMPORTANT NOTE ! !Note that convection is considered first, so that if !turbulence and convection occur in the same layer, the !erosion rate for turbulence is selected. ! if (eros_choice) then !Enhanced erosion in convective layers if (Mc(i,k,j) .gt. mc_thresh) tmp2s = eros_scale_c !Enhanced erosion in turbulent layers if (diff_t(i,k,j).gt.diff_thresh .or. & diff_t(i,k,min(j+1,KDIM)).gt.diff_thresh) & tmp2s = eros_scale_t end if !for erosion choice if (ql_upd(i,k) .gt. qmin .or. qi_upd(i,k) .gt. qmin) then D_eros(i,k)=qa_upd(i,k) * tmp2s * dtcloud * qs(i,k,j) * & (1.-U(i,k)) / (qi_upd(i,k) + ql_upd(i,k)) if (pfull(i,k,j) .gt. 400.e02) then D_eros(i,k)=D_eros(i,k)+efact*D_eros(i,k)*((pfull(i,k,kdim)- & pfull(i,k,j))/(pfull(i,k,kdim)-400.e02)) else D_eros(i,k)=D_eros(i,k)+efact*D_eros(i,k) endif else D_eros(i,k) = 0. endif enddo enddo ! ! The next step is to analytically integrate the saturated volume ! fraction equation. This follows the Tiedtke approach ! ! The qa equation is written in the form: ! ! (13) dqa/dt = (1.-qa) * A - qa * B ! ! Note that over the physics time step, A, B are assumed to be ! constants. ! ! Defining qa(t) = qa0 and qa(t+dtcloud) = qa1, the analytic ! solution of the above equation is: ! ! (14) qa1 = qaeq - (qaeq - qa0) * exp (-(A+B)*dtcloud) ! ! where qaeq is the equilibrium cloud fraction that is approached ! with an time scale of 1/(A+B), ! ! (15) qaeq = A/(A+B) ! ! ! To diagnose the magnitude of each of the right hand terms of ! (13) integrated over the time step, define the average cloud ! fraction in the interval t to t + dtcloud qabar as: ! ! (16) qabar = qaeq - [ (qa1-qa0) / ( dtcloud * (A+B) ) ] ! ! from which the magnitudes of the A and B terms integrated ! over the time step are: ! ! A * (1-qabar) and -B * (qabar) ! ! Additional notes on this analytic integration: ! ! 1. For large-scale cloud formation or destruction from ! the dqs_ls term the contributions to A or B are defined ! from: ! ! (19) A_ls * (1. - qa) = da_ls / dtcloud if da_ls >= 0. ! ! (20) B_ls * qa = da_ls / dtcloud if da_ls < 0. ! ! 2. Note that to reduce the number of variables, the following ! equivalency exists: ! ! Ql or Qi equation Qa equation ! -------------------- ------------------- ! ! C_dt A_dt ! D_dt B_dt ! qceq qaeq ! qcbar qabar ! qc1 qa1 ! qc0 qa0 ! ! 3. Qa goes to zero only in the case of ql and qi less than or ! equal to qmin; see 'cloud destruction code' near the end of ! this loop over levels. ! do k=1,jdim do i=1,idim !compute C_dt; This is assigned to the large-scale source term !following (18). Reset D_dt. C_dts = da_ls(i,k)/max((1.-qa_upd(i,k)),qmin) D_dts = D_eros(i,k) !do analytic integration qc0s = qa_upd(i,k) if ( (C_dts.gt.Dmin) .or. (D_dts.gt.Dmin) ) then qceqs = C_dts / (C_dts + D_dts) qc1s = qceqs - (qceqs - qc0s) * exp ( -1.*(C_dts+D_dts) ) qcbars = qceqs - ((qc1s - qc0s)/ (C_dts + D_dts)) else qceqs = qc0s qc1s = qc0s qcbars = qc0s endif !set total tendency term and update cloud fraction if (limit_conv_cloud_frac) then !RSH limit cloud area to be no more than that which is not being ! taken by convective clouds qc1s = MIN(qc1s, 1.0 -ahuco(i,k,j)) endif SA(i,k,j) = SA(i,k,j) + qc1s - qc0s qa_upd(i,k) = qc1s ! if (max(diag_id%qadt_lsform,diag_id%qa_lsform_col) > 0) qadt_lsform(i,k,j) = C_dts * (1.-qcbars) * inv_dtcloud if (max(diag_id%qadt_lsform,diag_id%qa_lsform_col) > 0) diag_4d(i,k,j,diag_pt%qadt_lsform) = C_dts * (1.-qcbars) * inv_dtcloud ! if (max(diag_id%qadt_eros,diag_id%qa_eros_col) > 0) qadt_eros (i,k,j) = D_dts * qcbars * inv_dtcloud if (max(diag_id%qadt_eros,diag_id%qa_eros_col) > 0) diag_4d(i,k,j,diag_pt%qadt_eros) = D_dts * qcbars * inv_dtcloud delta_cf(i,k) = C_dts * (1.-qcbars) ! The next step is to calculate the change in condensate ! due to non-convective condensation, dcond_ls. Note that this is ! not the final change but is used only to apportion condensate ! change between phases. According to Tiedtke 1993 this takes the ! form: ! ! (21) dcond_ls = -1. * (qa + 0.5*da_ls) * dqs_ls ! ! Here the 0.5*da_ls represents using a midpoint cloud fraction. ! This is accomplished by using the variable qcbar. dcond_ls(i,k) = -1. * qcbars * dqs_ls(i,k) enddo enddo !----------------------------------------------------------------------! ! ! ! ! ! ! ! NON-CONVECTIVE CONDENSATION ! ! ! ! ! ! ! ! METHOD 2 ! ! ! ! ! ! STATISTICAL CLOUD FRACTION ! ! ! ! ! else !for doing PDF cloud scheme !set Tiedtke erosion term to zero D_eros = 0. !compute pdf cloud fraction and condensate ! !Note that the SYMMETRIC beta distribution is used here. ! ! ! Initialize grid-box mean values of cloud fraction (qag), ! cloud condensate(qcg), and clear sky water vapor (qvg) qcg = 0. qvg = 0. qag = 0. !Create loop over sub-levels within a grid box do ns = 1, nsublevels !calculate normalized vertical level ! 0. = top of gridbox ! 1. = bottom of gridbox pnorm = (real(ns) - 0.5 )/real(nsublevels) !First step is to calculating the minimum (qtmin) !of the total water distribution and !the width of the qt distribution (deltaQ) ! !For diagnostic variance this is set to (1.-qthalfwidth)*qtbar !and 2*qthalfwidth*qtbar, respectively, where qtbar is the !mean total water in the grid box. ! ! qtbar = qta4(2,:,:,j)+pnorm*( (qta4(3,:,:,j)-qta4(2,:,:,j)) + & qta4(4,:,:,j)*(1-pnorm) ) qtbar = max(qmin,qtbar) deltaQ = 2.*qthalfwidth*qtbar qtmin = (1.-qthalfwidth)*qtbar !From this the variable normalized saturation specific !humidity qs_norm is calculated. ! ! qs_norm = (qs(Tl) - qtmin)/(qtmax-qtmin) ! ! = 0.5 - (qtbar - qs(Tl))/deltaQ ! !Note that if qs_norm > 1., the grid box is fully clear. !If qs_norm < 0., the grid box is fully cloudy. qs_norm = qtqsa4(2,:,:,j)+ & pnorm*( (qtqsa4(3,:,:,j)-qtqsa4(2,:,:,j)) + & qtqsa4(4,:,:,j)*(1-pnorm) ) qs_norm = 0.5 - ( qs_norm/deltaQ ) !Calculation of cloud fraction (qagtmp), cloud condensate !(qcgtmp), and water vapor in clear air part of the grid !box (qvgtmp) ! !Formulas (from Tompkins, and personal derivations): ! ! Define icbp = incomplete_beta(qs_norm,p,q) ! icbp1 = incomplete_beta(qs_norm,p+1,q) ! ! qagtmp = 1. - icbp ! ! qcgtmp = aThermo * { (qtbar-qtmin)*(1.-icbp1) - ! qs_norm*deltaQ*(1.-icbp ) } ! ! ! qvgtmp = qtmin + (p/(p+q))*(icbp1/icbp)*deltaQ ! ! ! where aThermo = 1./(1.+(L/cp)*dqsdT) ! ! note that in the qvg formula below the factor of 0.5 ! is equal to (p/(p+q)). ! do jd = 1,jdim do id = 1,idim if (qs_norm(id,jd).le.1.) then icbp = incomplete_beta(max(0.,qs_norm(id,jd)), & p = betaP , q = betaP) icbp1= incomplete_beta(max(0.,qs_norm(id,jd)), & p = betaP + 1, q = betaP) qagtmps = 1.-icbp qcgtmps = (qtbar(id,jd)-qtmin(id,jd))*(1.-icbp1)& - qs_norm(id,jd)*deltaQ(id,jd)*(1.-icbp) qcgtmps = qcgtmps/(1.+gamma(id,jd,j)) qvgtmps = qtmin(id,jd) + & 0.5*(icbp1/max(icbp,qmin))*deltaQ(id,jd) !bound very very small cloud fractions which may !cause negative cloud condensates due to roundoff !errors or similar errors in the beta table lookup. if((qagtmps.lt.0.).or.(qcgtmps.le.0.))then qagtmps = 0. qcgtmps = 0. qvgtmps = qtbar(id,jd) end if else qagtmps = 0. qcgtmps = 0. qvgtmps = qtbar(id,jd) end if !sum vertically ! !note special averaging of clear-sky water vapor !this is weighting clear-sky relative humidity by the !clear-sky fraction qag(id,jd) = qag(id,jd) + qagtmps qcg(id,jd) = qcg(id,jd) + qcgtmps qvg(id,jd) = qvg(id,jd)+(1.-qagtmps)*min(max(qvgtmps/max(qmin, & (qtbar(id,jd)+((qs_norm(id,jd)-0.5)*deltaQ(id,jd)))),0.),1.) !compute grid-box average cloud fraction, cloud condensate !and water vapor if (nsublevels.gt.1 .and. ns.eq.nsublevels) then qag(id,jd) = qag(id,jd) / real(nsublevels) qcg(id,jd) = qcg(id,jd) / real(nsublevels) !note special averaging of clear-sky water vapor if ((1.-qag(id,jd)).gt.qmin) then qvg(id,jd) =qvg(id,jd)/real(nsublevels)/& (1.-qag(id,jd)) qvg(id,jd) =qvg(id,jd) * qs(id,jd,j) else qvg(id,jd) = qs(id,jd,j) end if elseif (nsublevels.eq.1) then ! for nsublevels = 1, qag and qcg already hold their ! final values qvg(id,jd) = qvgtmps end if enddo enddo enddo !for number of sublevels loop !do adjustment of cloud fraction !rab qc0 = qa(:,:,j) !rab qc1 = qag !set total tendency term and update cloud fraction SA(:,:,j) = SA(:,:,j) + qag - qa(:,:,j) qa_upd = qag ! if (max(diag_id%qadt_lsform,diag_id%qa_lsform_col) > 0) qadt_lsform(:,:,j) = max(qag-qa(:,:,j),0.) * inv_dtcloud if (max(diag_id%qadt_lsform,diag_id%qa_lsform_col) > 0) diag_4d(:,:,j,diag_pt%qadt_lsform) = max(qag-qa(:,:,j),0.) * inv_dtcloud ! if (max(diag_id%qadt_lsdiss,diag_id%qa_lsdiss_col) > 0) qadt_lsdiss(:,:,j) = max(qa(:,:,j)-qag,0.) * inv_dtcloud if (max(diag_id%qadt_lsdiss,diag_id%qa_lsdiss_col) > 0) diag_4d(:,:,j,diag_pt%qadt_lsdiss) = max(qa(:,:,j)-qag,0.) * inv_dtcloud !define da_ls and delta_cf needed when do_liq_num = .true. (cjg) da_ls = max(qag-qa(:,:,j),0.) delta_cf = max(qag-qa(:,:,j),0.) !compute large-scale condensation / evaporation dcond_ls = qcg - (ql_upd + qi_upd) end if !for doing PDF clouds ! The next step is the apportionment on the non-convective ! condensation between liquid and ice phases. Following the ! suggestion of Rostayn (2000), all condensation at temperatures ! greater than -40C is in liquid form as ice nuclei are generally ! limited in the atmosphere. The droplets may subsequently be ! converted to ice by the Bergeron-Findeisan mechanism. ! ! One problem with this formulation is that the proper saturation ! vapor pressure is not used for cold clouds as it should be ! liquid saturation in the case of first forming liquid, but ! change to ice saturation as the cloud glaciates. The current ! use of ice saturation beneath -20C thus crudely mimics the ! result that nearly all stratiform clouds are glaciated for ! temperatures less than -15C. ! ! In the case of large-scale evaporation (dcond_ls<0.), it is ! assumed that cloud liquid will evaporate faster than cloud ! ice because if both are present in the same volume the ! saturation vapor pressure over the droplet is higher than ! that over the ice crystal. ! ! The fraction of large-scale condensation that is liquid ! is stored in the temporary variable tmp1. do k=1,jdim do i=1,idim !assume liquid fractionation tmp1s = 1. if (dcond_ls(i,k) .ge. 0.) then !For cases of cloud condensation where temperatures are !less than -40C create only ice if (T(i,k,j) .lt. tfreeze-40.) then tmp1s = 0. endif else if (qi_upd(i,k).gt.qmin) then if (ql_upd(i,k).gt.qmin) then !For cases of cloud evaporation of mixed phase clouds !set liquid evaporation to preferentially occur first tmp1s = min(-1.*dcond_ls(i,k),ql_upd(i,k))/max(-1.*dcond_ls(i,k),qmin) else !do evaporation of pure ice cloud tmp1s = 0. endif endif endif !calculate partitioning among liquid and ice to dcond_ls dcond_ls_ice(i,k) = (1.-tmp1s) * dcond_ls(i,k) dcond_ls(i,k) = tmp1s * dcond_ls(i,k) enddo enddo ! The next step is to compute semi-implicit qa,ql,qi which are ! used in many of the formulas below. This gives a somewhat ! implicitness to the scheme. In this calculation an estimate ! is made of what the cloud fields would be in the absence of ! cloud microphysics and cloud erosion. ! ! In the case of the Tiedtke cloud scheme, the mean cloud ! condensate is incremented if large-scale condensation is ! occuring. For cloud fraction, the value from the analytic ! integration above is used. ! ! For the statistical cloud scheme these are set equal to the ! values diagnosed from the beta-distribution apart from the ! corrections for mixed phase clouds. qa_mean = qa_upd if (.not. do_pdf_clouds) then ql_mean = ql_upd + max(dcond_ls ,0.) qi_mean = qi_upd + max(dcond_ls_ice,0.) else ql_mean = max(ql_upd + dcond_ls ,qmin) qi_mean = max(qi_upd + dcond_ls_ice,qmin) end if if (do_liq_num) then !yim's CCN activation do k = 1,jdim do i = 1,idim if ( da_ls(i,k) > 0.0 ) then up_strat = -1.*(((omega(i,k,j)+grav*Mc(i,k,j))/ & airdens(i,k,j)/grav) + & radturbten2(i,k,j)*cp_air/grav) !-->cjg: modification ! call aer_ccn_act (T(i,k,j), pfull(i,k,j), up_strat, & ! totalmass1(i,k,j,:), drop1(i,k)) thickness = deltpg(i,k) / airdens(i,k,j) wp2 = 2.0/(3.0*0.548**2)* & (0.5*(diff_t(i,k,j) + diff_t(i,k,min(j+1,KDIM)))/& thickness )**2 wp2 = MAX (wp2, var_limit**2) !rab take care of it when writing diags.... !rab debug2(i,k,j) = wp2**0.5 ! if (diag_id%debug2_3d > 0) debug2(i,k,j) = wp2 if (diag_id%subgrid_w_variance > 0 .and. up_strat > 0.0) diag_4d(i,k,j,diag_pt%subgrid_w_variance) = wp2 call aer_ccn_act_wpdf (T(i,k,j), pfull(i,k,j), & up_strat, wp2, & totalmass1(i,k,j,:), drop1(i,k)) ! if (diag_id%debug3_3d > 0) debug3(i,k,j) = drop1(i,k) if (diag_id%potential_droplets > 0 .and. up_strat > 0.0) diag_4d(i,k,j,diag_pt%potential_droplets) = drop1(i,k) ! if (diag_id%debug2_3d > 0) debug2(i,k,j) = 1. !<--cjg: end of modification qn_mean(i,k) = qn_upd(i,k) + max(delta_cf(i,k),0.)* & drop1(i,k)*1.e6/airdens(i,k,j) else drop1(i,k) = 0. qn_mean(i,k) = qn_upd(i,k) if (diag_id%subgrid_w_variance > 0) diag_4d(i,k,j,diag_pt%subgrid_w_variance) = 0.0 if (diag_id%potential_droplets > 0) diag_4d(i,k,j,diag_pt%potential_droplets) = 0.0 endif end do end do endif !compute diagnostics for cloud fraction ! if (diag_id%aall > 0) areaall(:,:,j) = qa_mean if (diag_id%aall > 0) diag_4d(:,:,j, diag_pt%aall) = qa_mean ! if (diag_id%aliq > 0 .or. diag_id%rvolume > 0) then if (max(diag_id%aliq,diag_id%rvolume) > 0) then ! where (ql_mean .gt. qmin) arealiq(:,:,j) = qa_mean where (ql_mean .gt. qmin) diag_4d(:,:,j,diag_pt%aliq) = qa_mean end if ! if (diag_id%aice > 0 .or. diag_id%vfall > 0) then if (max(diag_id%aice,diag_id%vfall) > 0) then ! where (qi_mean .gt. qmin) areaice(:,:,j) = qa_mean where (qi_mean .gt. qmin) diag_4d(:,:,j,diag_pt%aice) = qa_mean end if !----- -----! ! ! ! END OF ! ! ! ! NON-CONVECTIVE CONDENSATION ! ! ! ! SECTION ! ! ! ! ! ! ! !----------------------------------------------------------------------! !----------------------------------------------------------------------- ! ! ! ! TRANSFER OF RAIN FLUXES AND SNOW FLUXES BETWEEN ! ! SATURATED AND UNSATURATED PORTIONS OF THE GRID BOX ! ! AT LAYER INTERFACES ! ! ! Do transfer of rain and snow fluxes between clear and cloud ! portions of the grid box. This formulism follows the rain/snow ! parameterization developed by Jakob and Klein (2000). It ! treats the grid mean flux of rain and snow separately according ! to the portion that exists in unsaturated air and the portion ! that exists in saturated air. At the top of each level, some ! precipitation that was in unsaturated air in the levels above ! may enter saturated air and vice versa. These transfers are ! calculated here under an assumption for the overlap between ! precipitation area and saturated area at the same and adjacent ! levels. ! ! For rain, the area of the grid box in which rain in saturated ! air of the level above enters unsaturated air in the current ! level is called da_cld2clr and is given by ! maximum(a_rain_cld - qa , 0.) in the case of maximum overlap of ! rain_cld and qa, but equal to a_rain_cld*(1.-qa) in the case of ! random overlap of cloud and precipitation. The grid mean flux of ! rain which is transfered from saturated air to unsaturated air ! is given by rain_cld*da_cld2clr/a_rain_cld. ! ! The area of the grid box where rain in unsaturated air of the ! level above enters saturated air in the current level is ! called da_clr2cld and is given by ! maximum(0.,mininum(a_rain_clr,qa(j)-qa(j-1))) in the case of ! maximum overlap of cloud and rain_clr, but equal to a_rain_clr* ! qa for the case of random overlap of cloud and precipitation. ! The grid mean flux of rain transfered from unsaturated air to ! saturated air is given by rain_clr*da_clr2cld/a_rain_clr. ! ! NOTE: the overlap assumption used is set by the namelist ! variable in cloud_rad ! ! Overlap values are: 1 = maximum ! 2 = random ! ! If cloud_generator is on, the overlap choice is taken from ! there, by computing a quantity alpha, which is weighting ! between maximum and random overlap solutions. ! ! alpha = 1 ---> maximum, ! alpha = 0 ---> random ! ! alpha is stored in tmp3. ! tmp1 has the maximum overlap solution ! tmp2 has the random overlap solution ! if (cloud_generator_on) then if (j.gt.1) then tmp3 = compute_overlap_weighting(qa_mean_lst,qa_mean,& pfull(:,:,j-1),pfull(:,:,j)) tmp3 = min(1.,max(0.,tmp3)) else tmp3 = 0.0 end if end if ! ! ! Rain transfers are done first ! call cloud_clear_xfer (cloud_generator_on, tmp3, qa_mean, qa_mean_lst, a_rain_clr, a_rain_cld, rain_clr, rain_cld) ! ! Snow transfers are done second, in a manner exactly like that ! done for the rain fluxes ! call cloud_clear_xfer (cloud_generator_on, tmp3, qa_mean, qa_mean_lst, a_snow_clr, a_snow_cld, snow_clr, snow_cld) !----------------------------------------------------------------------- ! ! ! MELTING OF CLEAR SKY SNOW FLUX ! ! ! Melting of falling ice to rain occurs when T > tfreeze. The ! amount of melting is limited to the melted amount that would ! cool the temperature to tfreeze. ! ! In the snowmelt bug version, the temperature of melting was ! tfreeze + 2. like the original Tiedtke (1993) paper, instead of ! tfreeze. snow_fact=0. if (do_old_snowmelt) snow_fact=2. do k=1,jdim do i=1,idim !compute grid mean change in snow flux to cool the !grid box to tfreeze and store in temporary variable tmp1 tmp1s = cp_air*(T(i,k,j)-tfreeze-snow_fact)*deltpg(i,k)*inv_dtcloud/hlf ! If snow_clr > tmp1, then the amount of snow melted is ! limited to tmp1, otherwise melt snow_clr. The amount ! melted is stored in tmp2 tmp2s = max(min(snow_clr(i,k),tmp1s),0.) ST(i,k,j) = ST(i,k,j) - hlf*tmp2s*dtcloud/deltpg(i,k)/cp_air rain_clr(i,k) = rain_clr(i,k) + tmp2s !raise a_rain_clr to a_snow_clr IF AND only IF melting occurs !and a_rain_clr < a_snow_clr if (tmp2s .gt. 0. .and. a_snow_clr(i,k) .gt. qmin) & a_rain_clr(i,k) = max(a_rain_clr(i,k),a_snow_clr(i,k)) ! If all of the snow has melted, then zero out a_snow_clr if (snow_clr(i,k).lt.tmp1s .and. a_snow_clr(i,k).gt.qmin) then snow_clr(i,k) = 0. a_snow_clr(i,k) = 0. else snow_clr(i,k) = snow_clr(i,k) - tmp2s endif ! if (max(diag_id%snow_melt,diag_id%snow_melt_col) > 0) snow_melt(i,k,j) = tmp2s/deltpg(i,k) if (max(diag_id%snow_melt,diag_id%snow_melt_col) > 0) diag_4d(i,k,j,diag_pt%snow_melt) = tmp2s/deltpg(i,k) !----------------------------------------------------------------------- ! ! ! MELTING OF CLOUDY SKY SNOW FLUX ! ! ! Melting of falling ice to rain occurs when T > tfreeze. The ! amount of melting is limited to the melted amount that would ! cool the temperature to tfreeze. ! if (.not.do_old_snowmelt) then !compute grid mean change in snow flux to cool the !grid box to tfreeze and store in temporary variable tmp1 ! !note that tmp1 already has the value of this variable !from the clear-sky melt calculation, so one does not need !to repeat the calculation here. ! !However, note that clear-sky snow melt may have already !reduced the temperature of the grid box - this snow melt is in !variable tmp2 from lines above. Thus the amount that one !can melt is less. tmp1s = tmp1s - tmp2s ! If snow_cld > tmp1, then the amount of snow melted is ! limited to tmp1, otherwise melt snow_cld. The amount ! melted is stored in tmp2 tmp2s = max(min(snow_cld(i,k),tmp1s),0.) ST(i,k,j) = ST(i,k,j) - hlf*tmp2s*dtcloud/deltpg(i,k)/cp_air rain_cld(i,k) = rain_cld(i,k) + tmp2s !raise a_rain_cld to a_snow_cld IF AND only IF melting occurs !and a_rain_cld < a_snow_cld if (tmp2s .gt. 0. .and. a_snow_cld(i,k) .gt. qmin) & a_rain_cld(i,k) = max(a_rain_cld(i,k),a_snow_cld(i,k)) ! If all of the snow has melted, then zero out a_snow_cld if (snow_cld(i,k).lt.tmp1s .and. a_snow_cld(i,k).gt.qmin) then snow_cld(i,k) = 0. a_snow_cld(i,k) = 0. else snow_cld(i,k) = snow_cld(i,k) - tmp2s endif ! if (max(diag_id%snow_melt,diag_id%snow_melt_col) > 0) snow_melt(i,k,j) = snow_melt(i,k,j) + & ! tmp2s/deltpg(i,k) if (max(diag_id%snow_melt,diag_id%snow_melt_col) > 0) diag_4d(i,k,j,diag_pt%snow_melt) = diag_4d(i,k,j,diag_pt%snow_melt) + & tmp2s/deltpg(i,k) end if !for snowmelt bugfix enddo enddo !----------------------------------------------------------------------! ! ! ! COMPUTE SLOPE FACTOR FOR ICE MICROPHYSICS ! ! [The following microphysics follows that of Rotstayn (1997)] ! The number concentration of ice crystals of diameter D in the ! SIZE interval D to D+dD is assumed to be ! distributed as in a Marshall Palmer distribution : ! ! (22) N(D)dD = Nof * Exp( - lamda_f * D) ! ! The slope factor and intercept are not assumed to be constant, ! but the slope factor is ! ! (23) lamda_f = 1.6X10^(3.+0.023(tfreeze-T)) ! ! Integration of (22) over all particle sizes with a constant ! density of ice crystals , rho_ice, and assumed spherical shape ! yields a relationship between the intercept parameter and qi ! ! (24) Nof = airdens*qi_local*(lamda_f^4)/pi*rho_ice ! ! ! For the calculation of riming and sublimation of snow, ! lamda_f is needed, so it is calculated here. ! ! Also qi_mean is updated here with the flux of ice that falls ! into the cloudy portion of the grid box from above. This permits ! the Bergeron and riming process to know about the ice that falls ! into the grid box in the same time step. !Increment qi_mean by the ice flux entering the !the grid box. To convert ice_flux to units of condensate by !multiply by dtcloud and dividing by the mass per unit area !of the grid box. Implicit here is the assumption that the !ice entering the cloud will be spread instantaneously over !all of the cloudy area. qi_mean = qi_mean + snow_cld*dtcloud/deltpg !snow falling into cloud reduces the amount that !falls out of cloud: a loss of cloud ice from settling !is defined to be positive ! if (max(diag_id%qidt_fall,diag_id%qi_fall_col) > 0) qidt_fall(:,:,j)= -1.*snow_cld/deltpg if (max(diag_id%qidt_fall,diag_id%qi_fall_col) > 0) diag_4d(:,:,j,diag_pt%qidt_fall)= -1.*snow_cld/deltpg !compute lamda_f lamda_f = 1.6 * 10**(3.+0.023*(tfreeze-T(:,:,j))) !----------------------------------------------------------------------! ! ! ! ! ! ! ! LIQUID PHASE MICROPHYSICS ! ! ! ! AND ! ! ! ! ANALYTIC INTEGRATION OF QL EQUATION ! ! ! ! ! ! ! ! Accretion ! ! The parameterization of collection of cloud liquid drops by ! rain drops follows the parameterization of Rotstayn (1997). ! The parameterization starts with the continous-collection ! equation of drops by a rain drop of diameter D, falling with ! speed V(D). The fall speed of rain drops is taken from ! Gunn and Kinzer(1949): ! ! (25) V(D) = 141.4 (m**0.5,s-1) * sqrt(D) * (rho_ref/airdens)**0.5 ! ! where D is the radius of the rain drops in meters. Here ! rho_ref is a reference air density. This formula is generally ! good for 1 mm < D < 4 mm. ! ! The distribution of rain drops by SIZE follows a Marshall- ! Palmer distribution: ! ! (26) N(D) dD = Nor * Exp (- lamda *D) ! ! where N(D)dD is the number of rain drops with SIZE D in the ! interval D to D+dD, Nor is the intercept (assumed fixed at ! 8E+04 (1/m*m*m*m). lamda is the slope intercept parameter ! and with (21) it can be shown that lamda is a function of ! the rain rate. ! ! With these assumptions the local rate of accretion of cloud ! liquid reduces to: ! ! (27) dl/dt_local = - CB*Eco*((rain_rate_local/dens_h2o)**(7/9))* ! ! (rho_ref/airdens)**(1/9) * ql_local ! ! where CB is the accretion constant: ! ! CB = 65.772565 [(m)**(-7/9)] * [(s)**(-2/9)] ! ! AND Eco is the collection efficiency of a cloud droplet by a ! rain droplet. A good fit to the Table 8.2 of Rogers and Yau ! (1988) for rain drops between SIZE 1mm and 3mm is: ! ! (28) Eco = rad_liq**2 / (rad_liq**2 + 20.5 microns**2) . ! ! In generalizing to grid mean conditions (27) becomes: ! ! (29) dl/dt = - (arain_cld/qa_mean) * CB * Eco * ! ! [(rain_cld/a_rain_cld/dens_h2o)**(7/9)] * ql ! ! Note that the very weak dependence on air density is ! neglected at this point. ! ! The particle sizes are computed from the following equation ! ! (30) rad_liq = (3*airdens*ql/(4*pi*liq_dens*qa*N)^(1/3) ! ! ! For numerical treatment we write (25) as: ! ! dl/dt = - D_acc * l ! ! and if we do so: ! ! (31) D_acc = (arain_cld/qa_mean) * CB * Eco * ! ! [(rain_cld/a_rain_cld/dens_h2o)**(7/9)] ! ! In the work below, D_acc is added to D1_dt, the first processes ! contributing to the depletion of cloud liquid in the analytic ! integration. D1_dt represents the conversion of liquid to rain. if (.not. do_liq_num) then !compute rad_liq. The constant below is equal to !1.E+06 * (3/4*pi)^(1/3), where the 1E+06 is !the factor to convert meters to microns. rad_liq= 620350.49 *( (airdens(:,:,j)*ql_mean/ & max(qa_mean,qmin)/N/dens_h2o)**(1./3.)) !do not let very small cloud fractions contribution to !autoconversion or accretion where (qa_mean .le. qmin) rad_liq = 0. else !yim The 1st place droplet number is used rad_liq= 620350.49 *( (ql_mean/max(qn_mean,qmin)/dens_h2o)**(1./3.)) !do not let very small cloud fractions contribution to !autoconversion or accretion where (qa_mean .le. qmin .or. qn_upd .le.qmin) rad_liq = 0. endif ! if (diag_id%rvolume > 0) rvolume(:,:,j) = rad_liq*arealiq(:,:,j) ! if (diag_id%rvolume > 0) rvolume(:,:,j) = rad_liq*diag_4d(:,:,j,diag_pt%aliq) ! if (diag_id%rvolume > 0) diag_4d(:,:,j,diag_pt%rvolume) = rad_liq*arealiq(:,:,j) if (diag_id%rvolume > 0) diag_4d(:,:,j,diag_pt%rvolume) = rad_liq*diag_4d(:,:,j,diag_pt%aliq) !compute accretion D term D1_dt = dtcloud * 65.772565 * (a_rain_cld/max(qa_mean,qmin))* & ( rad_liq*rad_liq / (rad_liq*rad_liq+20.5) ) * & ((rain_cld/max(a_rain_cld,qmin)/dens_h2o)**(7./9.)) ! if (max(diag_id%qldt_accr,diag_id%ql_accr_col) > 0) qldt_accr(:,:,j) = D1_dt if (max(diag_id%qldt_accr,diag_id%ql_accr_col) > 0) diag_4d(:,:,j,diag_pt%qldt_accr) = D1_dt ! Autoconversion ! ! The autoconversion parameterization follow that of Manton ! and Cotton (1977). This formula has been used in Chen and ! Cotton (1987) and is used in the CSIRO GCM (Rotstayn 1997) ! and the LMD GCM (Boucher, Le Treut and Baker 1995). In this ! formulation the time rate of change of grid mean liquid is ! ! (32) dl/dt= -CA * qa * [(ql/qa)^(7/3)] * [(N*dens_h2o)^(-1/3)] ! ! * H(rad_liq - rthresh) ! ! where N is the number of cloud droplets per cubic metre, ! rthresh is a particle radius threshold needed to for autoconv- ! ersion to occur, H is the Heaviside function, and CA is ! a constant which is: ! ! (33) CA = 0.104 * grav * Ec * (airdens)^(4/3) / mu ! ! where grav is gravitational acceleration, Ec is the collection ! efficiency, airdens is the density of air, and mu is the ! dynamic viscosity of air. This constant is evaluated ! ignoring the temperature dependence of mu and with a fixed ! airdens of 1 kg/m3. ! ! With Ec = 0.55 (standard choice - see references) ! grav = 9.81 m/(s*s) ! airdens = 1.00 kg air/(m*m*m) ! mu = 1.717 E-05 kg condensate/(m*s) ! ! CA = 32681. [(kg air)^(4/3)]/kg liq/m/s ! ! ! For numerical treatment we write (32) as: ! ! dl/dt = - D_aut * l ! ! and if we do so: ! ! (34) D_aut = CA * [(N*dens_h2o)^(-1/3)] * [(ql/qa)^(4/3)] * & ! H(r-rthresh) ! ! In the work below, D_aut is temporarily stored in the variable ! tmp1 before being added to D1_dt. D1_dt represents the ! conversion of liquid to rain. ! ! Following Rotstayn, autoconversion is limited to the amount that ! would reduce the local liquid cloud condensate to the critical ! value at which autoconversion begins. This limiter is likely to ! be invoked frequently and is computed from ! ! (35) D_dt = log( (rad_liq/rthresh)**3. ) ! ! This limiter is stored in tmp2. ! ! ! ! ------------------------------------------------- ! ! Khairoutdinov and Kogan (2000) Autoconversion ! ! Reference: Khairoutdinov, M. and Y. Kogan, 2000: A new cloud ! physics parameterization in a large-eddy simulation ! model of marine stratocumulus. Mon. Wea. Rev., 128, ! 229-243. ! ! If the namelist parameter use_kk_auto = true, then the ! Khairoutdinov and Kogan (KK) autoconversion parameterization ! is used in place of the Manton and Cotton formula described ! above. ! ! In SI units this formula is: ! ! (32A) dl/dt= -CA * qa * [(ql/qa)^(2.47)] * [(N)^(-1.79)] ! ! ! where N is the number of cloud droplets per cubic metre ! and CA is a constant which is: ! ! (33A) CA = 7.4188E+13 (kg condensate/kg air)**(-1.47) ! (# drops/meters**3)**(1.79) ! seconds**(-1) ! ! For numerical treatment we write (32A) as: ! ! dl/dt = - D_aut * l ! ! and if we do so: ! ! (34A) D_aut = CA * [(N)^(-1.79)] * [(ql/qa)^(1.47)] ! if (do_liq_num) then if (use_kk_auto) then !*************************yim's version based on Khairoutdinov and Kogan (2000) !The second place N is used tmp1 = dtcloud * 1350. * & (1.e-6*max(qn_mean*airdens(:,:,j),max(qa_mean,qmin)*N_min))**(-1.79)* & (ql_mean)**(1.47)*max(qa_mean,qmin)**(0.32) ! tmp1 = dtcloud * 1350. * & ! (1.e-6*max(qn_mean,N_min)*airdens(:,:,j))**(-1.79)* & ! (ql_mean)**(1.47)*max(qa_mean,qmin)**(0.32) !************************** else !yim fall back to M & C using qn !compute autoconversion sink as in (34) tmp1 = 32681. * dtcloud * ((max(qn_mean,qmin)*airdens(:,:,j)*dens_h2o)**(-1./3.))* & (ql_mean**(4./3.))/max(qa_mean,qmin) !compute limiter as in (35) tmp2 =max(3*log(max(rad_liq,qmin)/rthresh),0.) !limit autoconversion to the limiter tmp1 = min(tmp1,tmp2) endif else if ( use_kk_auto ) then tmp1 = 0. !compute autoconversion sink as in (34A) where (ql_mean.gt.qmin) tmp1 = 7.4188E+13 * dtcloud * (N**(-1.79))* & ((ql_mean/max(qa_mean,qmin))**(1.47)) endwhere else !compute autoconversion sink as in (34) tmp1 = 32681. * dtcloud * ((N*dens_h2o)**(-1./3.))* & ((ql_mean/max(qa_mean,qmin))**(4./3.)) !compute limiter as in (35) tmp2 =max(3*log(max(rad_liq,qmin)/rthresh),0.) !limit autoconversion to the limiter tmp1 = min(tmp1,tmp2) endif endif !add autoconversion to D1_dt D1_dt = D1_dt + tmp1 !auto conversion will change a_rain_cld upto area of cloud where (tmp1 .gt. Dmin) a_rain_cld = qa_mean ! if (max(diag_id%qldt_auto,diag_id%ql_auto_col) > 0) qldt_auto(:,:,j) = tmp1 if (max(diag_id%qldt_auto,diag_id%ql_auto_col) > 0) diag_4d(:,:,j,diag_pt%qldt_auto) = tmp1 if ( diag_id%aauto > 0 ) then ! where ( rad_liq .gt. rthresh ) areaautocv(:,:,j) = qa_mean where ( rad_liq .gt. rthresh ) diag_4d(:,:,j,diag_pt%aauto) = qa_mean end if ! Bergeron-Findeisan Process ! ! where ice and liquid coexist, the differential saturation ! vapor pressure between liquid and ice phases encourages ! the growth of ice crystals at the expense of liquid droplets. ! ! Rotstayn (2000) derive an equation for the growth of ice by ! starting with the vapor deposition equation for an ice crystal ! and write it in terms of ice specific humidity as: ! ! {(Ni/airdens)**2/3}*7.8 ! (36) dqi/dt = ------------------------ X [(esl-esi)/esi] X ! [rhoice**1/3]* A_plus_B ! ! ((max(qi,Mio*Ni/airdens))**(1/3))* ! ! Here Ni is the ice crystal number which is taken from the ! parameterization of Meyers et al. : ! ! (37) Ni = 1000 * exp( (12.96* [(esl-esi)/esi]) - 0.639 ) ! ! The use of the maximum operator assumed that there is a ! background ice crystal always present on which deposition can ! occur. Mio is an initial ice crystal mass taken to be 10-12kg. ! ! Figure 9.3 of Rogers and Yau (1998) shows the nearly linear ! variation of [(esl-esi)/esi] from 0. at 273.16K to 0.5 at ! 233.16K. Analytically this is parameterized as (tfreeze-T)/80. ! ! ! Generalizing (36) to grid mean conditions and writing it in ! terms of a loss of cloud liquid yields: ! ! (1000*exp((12.96*(tfreeze-T)/80)-0.639)/airdens)**2/3 ! (38) dql/dt = - ----------------------------------------------------- ! [rhoice**1/3]* A_plus_B * ql ! ! *qa*7.8*((max(qi/qa,Mio*Ni/airdens))**(1/3))*[(tfreeze-T)/80]*ql ! ! Note that the density of ice is set to 700 kg m3 the value ! appropriate for pristine ice crystals. This value is ! necessarily different than the value of ice used in the riming ! and sublimation part of the code which is for larger particles ! which have much lower densities. if (do_dust_berg) then crystal=0. do k = 1,jdim do i = 1,idim if ( (T(i,k,j) .lt. tfreeze) .and. (ql_mean(i,k) .gt. qmin) & .and. (qa_mean(i,k) .gt. qmin)) then Si0=1+0.0125*(tfreeze-T(i,k,j)) call Jhete_dep(T(i,k,j),Si0,concen_dust_sub(i,k,j),crystal(i,k)) ! if (diag_id%debug4_3d > 0) debug4(i,k,j) = 1. if (diag_id%dust_berg_flag > 0) diag_4d(i,k,j,diag_pt%dust_berg_flag) = 1. endif end do end do ! if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) qndt_cond(:,:,j) = crystal if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) diag_4d(:,:,j, diag_pt%qndt_cond) = crystal if (max(diag_id%qndt_evap,diag_id%qn_evap_col) > 0) then ! qndt_evap(:,:,j) = 0. diag_4d(:,:,j,diag_pt%qndt_evap) = 0. where (T(:,:,j).lt.tfreeze .and. ql_mean.gt.qmin .and. qa_mean.gt.qmin) ! qndt_evap(:,:,j) = 1.e-3*exp((12.96*0.0125*(tfreeze-T(:,:,j)))-0.639) diag_4d(:,:,j,diag_pt%qndt_evap) = 1.e-3*exp((12.96*0.0125*(tfreeze-T(:,:,j)))-0.639) end where endif !do Bergeron process D2_dt = 0.0 where (T(:,:,j) .lt. tfreeze .and. ql_mean .gt. qmin .and. qa_mean .gt. qmin) D2_dt = dtcloud * qa_mean * ((1.e6*crystal(:,:)/airdens(:,:,j))**(2./ & 3.))* 7.8* ((max(qi_mean/qa_mean,1.E-12*1.e6* & crystal(:,:) & /airdens(:,:,j)))**(1./3.))*0.0125* & (tfreeze-T(:,:,j))/((700.**(1./3.))* & A_plus_B(:,:,j)*ql_mean) end where else !do Bergeron process D2_dt = 0.0 where (T(:,:,j) .lt. tfreeze .and. ql_mean .gt. qmin .and. qa_mean .gt. qmin) D2_dt = dtcloud * qa_mean * ((cfact*1000.*exp((12.96*0.0125* & (tfreeze-T(:,:,j)))-0.639)/airdens(:,:,j))**(2./ & 3.))* 7.8* ((max(qi_mean/qa_mean,1.E-12*cfact*1000.* & exp((12.96*0.0125*(tfreeze-T(:,:,j)))-0.639) & /airdens(:,:,j)))**(1./3.))*0.0125* & (tfreeze-T(:,:,j))/((700.**(1./3.))* & A_plus_B(:,:,j)*ql_mean) end where endif sum_berg = D2_dt ! if (max(diag_id%qldt_berg,diag_id%ql_berg_col) > 0) qldt_berg(:,:,j) = D2_dt if (max(diag_id%qldt_berg,diag_id%ql_berg_col) > 0) diag_4d(:,:,j,diag_pt%qldt_berg) = D2_dt ! Accretion of cloud liquid by ice ('Riming') ! ! [The below follows Rotstayn] ! Accretion of cloud liquid by ice ('Riming') is derived in ! the same way as the accretion of cloud liquid by rain. That ! is the continous-collection equation for the growth of an ! ice crystal is integrated over all ice crystal sizes. This ! calculation assumes all crystals fall at the mass weighted ! fall speed, Vfall. This yields the following equation after ! accounting for the area of the interaction ! ! (39) dql/dt = - ( a_snow_cld / qa/a_snow_cld ) * ! ( ELI * lamda_f * snow_cld / 2/ rho_ice ) * ql ! ! ! Note that in the version with the snowmelt bug, riming was ! prevented when temperatures were in excess of freezing. !add in accretion of cloud liquid by ice tmp1 = 0.0 if (do_old_snowmelt) then where ((a_snow_cld.gt.qmin) .and. (ql_mean.gt.qmin) .and. & ( qa_mean.gt.qmin) .and. (T(:,:,j) .lt. tfreeze) ) tmp1 = dtcloud*0.5*ELI*lamda_f*snow_cld/qa_mean/rho_ice end where else where ((a_snow_cld.gt.qmin) .and. (ql_mean.gt.qmin) .and. & ( qa_mean.gt.qmin) ) tmp1 = dtcloud*0.5*ELI*lamda_f*snow_cld/qa_mean/rho_ice end where end if D2_dt = D2_dt + tmp1 sum_rime = tmp1 ! if (max(diag_id%qldt_rime,diag_id%ql_rime_col) > 0) qldt_rime(:,:,j) = tmp1 if (max(diag_id%qldt_rime,diag_id%ql_rime_col) > 0) diag_4d(:,:,j,diag_pt%qldt_rime) = tmp1 ! Freezing of cloud liquid to cloud ice occurs when ! the temperature is less than -40C. At these very cold temper- ! atures it is assumed that homogenous freezing of cloud liquid ! droplets will occur. To accomplish this numerically in one ! time step: ! ! (40) D*dtcloud = ln( ql / qmin ). ! ! With this form it is guaranteed that if this is the only ! process acting that ql = qmin after one integration. ! !do homogeneous freezing where ( T(:,:,j).lt.tfreeze-40..and.(ql_mean.gt.qmin).and. & (qa_mean.gt.qmin)) D2_dt = log ( ql_mean / qmin ) end where do k=1,jdim do i=1,idim if (T(i,k,j).lt.(tfreeze-40.).and.(ql_mean(i,k).gt.qmin) & .and.(qa_mean(i,k).gt.qmin)) then sum_freeze(i,k) = D2_dt(i,k) sum_rime(i,k) = 0. sum_berg(i,k) = 0. endif enddo enddo if (max(diag_id%qldt_freez,diag_id%qldt_rime,diag_id%qldt_berg,diag_id%ql_freez_col, & diag_id%ql_rime_col,diag_id%ql_berg_col) > 0) then do k=1,jdim do i=1,idim if (T(i,k,j).lt.(tfreeze-40.).and.(ql_mean(i,k).gt.qmin) & .and.(qa_mean(i,k).gt.qmin)) then if (max(diag_id%qldt_freez,diag_id%ql_freez_col) > 0) diag_4d(i,k,j,diag_pt%qldt_freez) = D2_dt(i,k) if (max(diag_id%qldt_rime,diag_id%ql_rime_col) > 0) diag_4d(i,k,j,diag_pt%qldt_rime ) = 0. if (max(diag_id%qldt_berg,diag_id%ql_berg_col) > 0) diag_4d(i,k,j,diag_pt%qldt_berg ) = 0. endif enddo enddo end if ! For BC aerosol in-cloud scavenging: do k=1,jdim do i=1,idim qldt_sum = sum_berg(i,k) + sum_freeze(i,k) + sum_rime(i,k) if (qldt_sum > 0.) f_snow_berg(i,k,j) = sum_berg(i,k)/qldt_sum enddo enddo ! Analytic integration of ql equation ! ! ! The next step is to analytically integrate the cloud liquid ! condensate equation. This follows the Tiedtke approach. ! ! The qc equation is written in the form: ! ! (41) dqc/dt = C - qc * D ! ! Note that over the physics time step, C and D are assumed to ! be constants. ! ! Defining qc(t) = qc0 and qc(t+dtcloud) = qc1, the analytic ! solution of the above equation is: ! ! (42) qc1 = qceq - (qceq - qc0) * exp (-D*dtcloud) ! ! where qceq is the equilibrium cloud condensate that is approached ! with an time scale of 1/D, ! ! (43) qceq = C / D ! ! ! To diagnose the magnitude of each of the right hand terms of ! (41) integrated over the time step, define the average cloud ! condensate in the interval t to t + dtcloud qcbar as: ! ! (44) qcbar = qceq - [ (qc1-qc0) / ( dtcloud * D ) ] ! ! from which the magnitudes of the C and D terms integrated ! over the time step are: ! ! C and -D * (qcbar) ! ! ! Additional notes on this analytic integration: ! ! 1. Because of finite machine precision it is required that ! D*dt is greater than a minimum value. This minimum ! alue occurs where 1. - exp(-D*dt) = 0. instead of ! D*dt. This value will be machine dependent. See discussion ! at top of code for Dmin. ! !C_dt is set to large-scale condensation. Sink of cloud liquid !is set to the sum of D1 (liquid to rain component), and D2 !(liquid to ice component), D_eros (erosion), and large-scale !evaporation (note use of ql mean). do k=1,jdim do i=1,idim C_dts = max(dcond_ls(i,k),0.) D_dts = D1_dt(i,k) + D2_dt(i,k) + D_eros(i,k) + & (max(-1.*dcond_ls(i,k),0.)/max(ql_mean(i,k),qmin)) !do analytic integration qc0s = ql_upd(i,k) if ( D_dts.gt.Dmin ) then qceqs = C_dts / D_dts qc1s = qceqs - (qceqs - qc0s) * exp ( -1.* D_dts ) qcbars = qceqs - ((qc1s - qc0s)/ D_dts) else qceqs = qc0s + C_dts qc1s = qc0s + C_dts qcbars = qc0s + 0.5*C_dts endif !set total tendency term and update cloud !Note that the amount of SL calculated here is stored in tmp1. SL(i,k,j) = SL(i,k,j) + qc1s - qc0s ql_upd(i,k) = qc1s !compute the amount each term contributes to the change !rab Dterm = -D_dt * qcbar ! Apportion SL between various processes. This is necessary to ! account for how much the temperature changes due to various ! phase changes. For example: ! ! liquid to ice = (D2/D)*(-Dterm) ! ! liquid to rain = (D1/D)*(-Dterm) ! ! (no phase change but needed to know how much to increment ! rainflux) ! ! vapor to liquid = - { ((-dcond_ls/ql_mean)+D_eros)/D}*(-Dterm) ! where dcond_ls < 0 ! ! but ! ! dcond_ls -(D_eros/D)*(-Dterm) ! where dcond_ls > 0 ! !initialize tmp2 to hold (-Dterm)/D !rab tmp2 = -Dterm/max(D_dt,Dmin) if (Nml%retain_cm3_bug) then tmp2(i,k) = D_dts*qcbars/max(D_dts,Dmin) else if ( D_dts.gt.Dmin ) then tmp2(i,k) = D_dts*qcbars/max(D_dts,Dmin) else tmp2(i,k) = 0. endif endif !do phase changes from large-scale processes and boundary !layer condensation/evaporation ST(i,k,j) = ST(i,k,j) + (hlv*max(dcond_ls(i,k),0.)/cp_air) - & (hlv*(max(-1.*dcond_ls(i,k),0.) /max(ql_mean(i,k),qmin))*tmp2(i,k)/cp_air) SQ(i,k,j) = SQ(i,k,j) - max(dcond_ls(i,k),0.) + & (max(-1.*dcond_ls(i,k),0.) /max(ql_mean(i,k),qmin))*tmp2(i,k) !add in liquid to ice and cloud erosion to temperature tendency ST(i,k,j) = ST(i,k,j) + (hlf*D2_dt(i,k)-hlv*D_eros(i,k))*tmp2(i,k)/cp_air !cloud evaporation adds to water vapor SQ(i,k,j) = SQ(i,k,j) + D_eros(i,k)*tmp2(i,k) !add conversion of liquid to rain to the rainflux rain_cld(i,k) = rain_cld(i,k) +D1_dt(i,k)*tmp2(i,k)*deltpg(i,k)*inv_dtcloud !save liquid converted to ice into tmp3 and increment qi_mean tmp3(i,k) = tmp2(i,k)*D2_dt(i,k) qi_mean(i,k) = qi_mean(i,k) + tmp3(i,k) enddo enddo if (do_liq_num) then !****************************************************************** ! Analytic integration of qn equation ! ! The qn equation is written in the form: ! ! (m1) dqc/dt = (1 - qabar) * A * qc^ - qc * D ! = C - qc * D ! ! where qc^ is the large-scale qn calculated from the ! activation parameterization. Note that over the physics ! time step, C and D are assumed to be constants. ! ! Defining qc(t) = qc0 and qc(t+dtcloud) = qc1, the analytic ! solution of the above equation is: ! ! (m2) qc1 = qceq - (qceq - qc0) * exp (-D*dtcloud) ! ! where qceq is the equilibrium cloud droplet number that is approached ! with an time scale of 1/D, ! ! (m3) qceq = C / D ! ! ! To diagnose the magnitude of each of the right hand terms of ! (m1) integrated over the time step, define the average cloud ! condensate in the interval t to t + dtcloud qcbar as: ! ! (m4) qcbar = qceq - [ (qc1-qc0) / ( dtcloud * D ) ] ! ! from which the magnitudes of the C and D terms integrated ! over the time step are: ! ! C and -D * (qcbar) ! !Calculate C_dt ! C_dt=max(tmp5,0.)*drop1*1.e6/airdens(:,:,j) !For replying the review, substract autoconversion ! D_dt = D1_dt + D2_dt + D_eros !do analytic integration ! where ( (D_dt.gt.Dmin) ) ! qc0 = qn_upd ! qceq = C_dt / max(D_dt, Dmin) ! qc1 = qceq - (qceq - qc0) * exp ( -1.* D_dt ) ! qcbar = qceq - ((qc1 - qc0)/ max(D_dt, Dmin)) ! elsewhere ! qc0 = qn_upd ! qceq = qc0 + C_dt ! qc1 = qc0 + C_dt ! qcbar = qc0 + 0.5*C_dt ! end where ! if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) qndt_cond(:,:,j) = 0. if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) diag_4d(:,:,j,diag_pt%qndt_cond) = 0. do k=1,jdim do i=1,idim !Calculate C_dt C_dts=max(delta_cf(i,k),0.)*drop1(i,k)*1.e6/airdens(i,k,j) D_dts = num_mass_ratio1*D1_dt(i,k) + (num_mass_ratio2*D2_dt(i,k) + D_eros(i,k)) qc0s = qn_upd(i,k) if (D_dts > Dmin) then qceqs = C_dts / D_dts qc1s = qceqs - (qceqs - qc0s)* exp(-1.*D_dts) qcbars = qceqs - ((qc1s -qc0s)/D_dts) else qceqs = qc0s + C_dts qc1s = qc0s + C_dts qcbars = qc0s + 0.5*C_dts endif !set total tendency term and update cloud !Note that the amount of SN calculated here is stored in tmp1. SN(i,k,j) = SN(i,k,j) + qc1s - qc0s qn_upd(i,k) = qc1s !compute the amount each term contributes to the change ! where ( C_dt .gt. 0 ) ! Cterm = C_dt ! Dterm = D_dt * qcbar ! elsewhere ! Cterm = 0. ! Dterm = D_dt * qcbar ! end where if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) then ! if ( C_dts.gt. 0 ) qndt_cond(i,k,j) = C_dts if ( C_dts.gt. 0 ) diag_4d(i,k,j,diag_pt%qndt_cond) = C_dts endif ! if (max(diag_id%qndt_evap,diag_id%qn_evap_col) > 0) qndt_evap(i,k,j) = D_dts * qcbars !Dterm if (max(diag_id%qndt_evap,diag_id%qn_evap_col) > 0) diag_4d(i,k,j,diag_pt%qndt_evap) = D_dts * qcbars !Dterm end do end do endif ! (do_liq_num) !**************************************************************************** ! ! diagnostics for cloud liquid tendencies ! ! if (max(diag_id%qldt_cond,diag_id%ql_cond_col) > 0) qldt_cond(:,:,j) = max(dcond_ls,0.) *inv_dtcloud if (max(diag_id%qldt_cond,diag_id%ql_cond_col) > 0) diag_4d(:,:,j,diag_pt%qldt_cond) = max(dcond_ls,0.) *inv_dtcloud ! if (max(diag_id%qldt_evap,diag_id%ql_evap_col) > 0) qldt_evap(:,:,j) = (max(0.,-1.*dcond_ls )/max(ql_mean, & if (max(diag_id%qldt_evap,diag_id%ql_evap_col) > 0) diag_4d(:,:,j,diag_pt%qldt_evap) = (max(0.,-1.*dcond_ls )/max(ql_mean, & qmin)) *tmp2*inv_dtcloud ! if (max(diag_id%qldt_accr,diag_id%ql_accr_col) > 0) qldt_accr(:,:,j) = qldt_accr (:,:,j)*tmp2*inv_dtcloud if (max(diag_id%qldt_accr,diag_id%ql_accr_col) > 0) diag_4d(:,:,j,diag_pt%qldt_accr) = diag_4d(:,:,j,diag_pt%qldt_accr)*tmp2*inv_dtcloud ! if (max(diag_id%qldt_auto,diag_id%ql_auto_col) > 0) qldt_auto(:,:,j) = qldt_auto (:,:,j)*tmp2*inv_dtcloud if (max(diag_id%qldt_auto,diag_id%ql_auto_col) > 0) diag_4d(:,:,j,diag_pt%qldt_auto) = diag_4d(:,:,j,diag_pt%qldt_auto )*tmp2*inv_dtcloud ! if (max(diag_id%qldt_eros,diag_id%ql_eros_col) > 0) qldt_eros(:,:,j) = D_eros *tmp2*inv_dtcloud ! if (max(diag_id%qldt_berg,diag_id%ql_berg_col) > 0) qldt_berg(:,:,j) = qldt_berg (:,:,j)*tmp2*inv_dtcloud if (max(diag_id%qldt_berg,diag_id%ql_berg_col) > 0) diag_4d(:,:,j,diag_pt%qldt_berg) = diag_4d(:,:,j,diag_pt%qldt_berg )*tmp2*inv_dtcloud ! if (max(diag_id%qldt_rime,diag_id%ql_rime_col) > 0) qldt_rime(:,:,j) = qldt_rime (:,:,j)*tmp2*inv_dtcloud if (max(diag_id%qldt_rime,diag_id%ql_rime_col) > 0) diag_4d(:,:,j,diag_pt%qldt_rime) = diag_4d(:,:,j,diag_pt%qldt_rime )*tmp2*inv_dtcloud ! if (max(diag_id%qldt_freez,diag_id%ql_freez_col) > 0) qldt_freez(:,:,j) = qldt_freez(:,:,j)*tmp2*inv_dtcloud if (max(diag_id%qldt_freez,diag_id%ql_freez_col) > 0) diag_4d(:,:,j,diag_pt%qldt_freez) = diag_4d(:,:,j,diag_pt%qldt_freez)*tmp2*inv_dtcloud ! if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) qndt_cond(:,:,j) = qndt_cond(:,:,j)*inv_dtcloud if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) diag_4d(:,:,j,diag_pt%qndt_cond) = diag_4d(:,:,j,diag_pt%qndt_cond)*inv_dtcloud ! if (max(diag_id%qndt_evap,diag_id%qn_evap_col) > 0) qndt_evap(:,:,j) = qndt_evap(:,:,j)*inv_dtcloud if (max(diag_id%qndt_evap,diag_id%qn_evap_col) > 0) diag_4d(:,:,j,diag_pt%qndt_evap) = diag_4d(:,:,j,diag_pt%qndt_evap)*inv_dtcloud ! if (max(diag_id%qldt_cond,diag_id%ql_cond_col) > 0) diag_4d(:,:,j,diag_pt%qldt_cond) = max(dcond_ls,0.) *inv_dtcloud ! if (max(diag_id%qldt_evap,diag_id%ql_evap_col) > 0) diag_4d(:,:,j, diag_pt%qldt_evap) = (max(0.,-1.*dcond_ls )/max(ql_mean, & ! qmin)) *tmp2*inv_dtcloud ! if (max(diag_id%qldt_accr,diag_id%ql_accr_col) > 0) diag_4d(:,:,j, diag_pt%qldt_accr) = qldt_accr (:,:,j) ! if (max(diag_id%qldt_auto,diag_id%ql_auto_col) > 0) diag_4d(:,:,j, diag_pt%qldt_auto) = qldt_auto (:,:,j) if (max(diag_id%qldt_eros,diag_id%ql_eros_col) > 0) diag_4d(:,:,j, diag_pt%qldt_eros) = D_eros *tmp2*inv_dtcloud ! if (max(diag_id%qldt_berg,diag_id%ql_berg_col) > 0) diag_4d(:,:,j, diag_pt%qldt_berg) = qldt_berg (:,:,j) ! if (max(diag_id%qldt_rime,diag_id%ql_rime_col) > 0) diag_4d(:,:,j, diag_pt%qldt_rime) = qldt_rime (:,:,j) ! if (max(diag_id%qldt_freez,diag_id%ql_freez_col) > 0) diag_4d(:,:,j,diag_pt%qldt_freez) = qldt_freez(:,:,j) ! if (max(diag_id%qndt_cond,diag_id%qn_cond_col) > 0) diag_4d(:,:,j, diag_pt%qndt_cond) = qndt_cond(:,:,j) ! if (max(diag_id%qndt_evap,diag_id%qn_evap_col) > 0) diag_4d(:,:,j, diag_pt%qndt_evap) = qndt_evap(:,:,j) !----- -----! ! ! ! END OF ! ! ! ! LIQUID PHASE MICROPHYSICS ! ! ! ! AND ! ! ! ! ANALYTIC INTEGRATION OF QL EQUATION ! ! ! ! SECTION ! ! ! ! ! ! ! !----------------------------------------------------------------------! !----------------------------------------------------------------------! ! ! ! ! ! ! ! ICE PHASE MICROPHYSICS ! ! ! ! AND ! ! ! ! ANALYTIC INTEGRATION OF QI EQUATION ! ! ! ! ! ! ! ! Ice settling ! ! Ice settling is treated as in Heymsfield Donner 1990. ! The mass weighted fall speed is parameterized as in equation ! #49. ! ! In terms of the analytic integration with respect qi of Tiedtke, ! the flux in from the top of the grid layer is equated to the ! source term, and the flux out of the bottom of the layer is ! equated to the sink term: ! ! (47) C_dt = snow_cld * dtcloud * grav / deltp ! ! (48) D_dt = airdens * grav * Vfall * dtcloud / deltp ! ! All ice crystals are assumed to fall with the same fall speed ! which is given as in Heymsfield and Donner (1990) as: ! ! (49) Vfall = 3.29 * ( (airdens*qi_mean/qa_mean)**0.16) ! ! which is the formula in Heymsfield and Donner. Note however ! that because this is uncertain, sensitivity runs will be made ! with different formulations. Note that when Vfall is computed ! the source incremented qi, qi_mean, is used. This gives some ! implicitness to the scheme. !compute Vfall iwc = airdens(:,:,j)*qi_mean/max(qa_mean,qmin) where (iwc >= iwc_crit) Vfall = vfact*3.29 * iwc**0.16 elsewhere Vfall = vfact*vfall_const2 * iwc**vfall_exp2 end where ! if (diag_id%vfall > 0) vfalldiag(:,:,j) = Vfall(:,:)*areaice(:,:,j) ! if (diag_id%vfall > 0) vfalldiag(:,:,j) = Vfall(:,:)*diag_4d(:,:,j,diag_pt%aice) ! if (diag_id%vfall > 0) diag_4d(:,:,j,diag_pt%vfall) = Vfall(:,:)*areaice(:,:,j) if (diag_id%vfall > 0) diag_4d(:,:,j,diag_pt%vfall) = Vfall(:,:)*diag_4d(:,:,j,diag_pt%aice) !add to ice source the settling ice flux from above !also note that tmp3 contains the source !of liquid converted to ice from above tmp3 = tmp3 + snow_cld*dtcloud/deltpg !Compute settling of ice. The result is multiplied by !dtcloud/deltp to convert to units of D_dt. !Note that if tracers are not advected then this is done !relative to the local vertical motion. if (tracer_advec) then tmp1 = 0. else tmp1 = omega(:,:,j) end if where (qi_mean .gt. qmin .and. qa_mean .gt. qmin) D1_dt = max(0.,((airdens(:,:,j)*Vfall)+(tmp1/grav))* & dtcloud/deltpg ) a_snow_cld = qa_mean elsewhere D1_dt = 0. snow_cld = 0. a_snow_cld = 0. end where ! Melting of in-cloud ice ! ! Melting occurs where the temperature is greater than Tfreezing. ! This is an instaneous process such that no stratiform ice will ! remain in the grid at the end of the timestep. ! No ice settles out of the grid box (i.e. D1 is set to zero) when ! the amount of ice to be melted is less than that that would ! bring the grid box to the freezing point. ! ! The ice that melts becomes rain. This is because if an ice ! crystal of dimension ~100 microns and mass density of 100 kg/m2 ! melts it will become a droplet of SIZE 40 microns which is ! clearly a drizzle SIZE drop. Ice crystals at temperatures near ! freezing are assumed to be this large, consistent with the ! assumption of particle SIZE temperature dependence. ! !compute grid mean change in cloud ice to cool the !grid box to 0C and store in temporary variable tmp1 tmp1 = cp_air*(T(:,:,j)-tfreeze)/hlf ! If qi_mean > tmp1, then the amount of ice melted is ! limited to tmp1, otherwise melt all qi_mean. The amount ! melted is stored in tmp2 tmp2 = max(min(qi_mean,tmp1),0.) D2_dt = max(0.,log(max(qi_mean,qmin)/max(qi_mean-tmp2,qmin))) !melting of ice creates area to a_rain_cld where (D2_dt .gt. Dmin) a_rain_cld = qa_mean endwhere !If all of the ice can melt, then don't permit any ice to fall !out of the grid box and set a_snow_cld to zero. where (qi_mean .lt. tmp1 .and. qi_mean .gt. qmin) D1_dt = 0. snow_cld = 0. a_snow_cld = 0. end where ! Analytic integration of qi equation ! ! This repeats the procedure done for the ql equation. ! See above notes for detail. !At this point C_dt already includes the source of cloud ice !falling from above as well as liquid converted to ice. !Therefore add in large_scale deposition. ! !Compute D_dt which has contributions from D1_dt (ice settling) !and D2_dt (ice melting), D_eros (cloud erosion), and large- !scale sublimation (note use of qi mean). do k=1,jdim do i=1,idim C_dts = tmp3(i,k) + max(dcond_ls_ice(i,k),0.) D_dts = D1_dt(i,k) + D2_dt(i,k) + D_eros(i,k) + & (max(-1.*dcond_ls_ice(i,k),0.)/max(qi_mean(i,k),qmin)) !do analytic integration qc0s = qi_upd(i,k) if ( D_dts.gt.Dmin ) then qceqs = C_dts / D_dts qc1s = qceqs - (qceqs - qc0s) * exp ( -1.* D_dts ) qcbars = qceqs - ((qc1s - qc0s)/D_dts) else qceqs = qc0s + C_dts qc1s = qc0s + C_dts qcbars = qc0s + 0.5*C_dts endif !set total tendency term and update cloud !Note that the amount of SL calculated here is stored in tmp1. SI(i,k,j) = SI(i,k,j) + qc1s - qc0s qi_upd(i,k) = qc1s !compute the amount each term contributes to the change !rab Dterm = -D_dt * qcbar ! Apportion SI between various processes. This is necessary to ! account for how much the temperature and water vapor changes ! due to various phase changes. For example: ! ! ice settling = (D1/D)*(-Dterm)*deltp/grav/dtcloud ! ! vapor to ice = ! -{ ((-dcond_ls_ice/qi_mean)+D_eros)/ D }*(-Dterm) ! where dcond_ls_ice < 0. ! ! but ! ! dcond_ls_ice -(D_eros/D)* (-Dterm) ! ! where dcond_ls_ice > 0. ! ! melting of ice = (D2/D)*(-Dterm)*deltp/grav/dtcloud ! !initialize tmp2 to hold (-Dterm)/D !rab tmp2 = -Dterm/max(D_dt,Dmin) if (Nml%retain_cm3_bug) then tmp2s = D_dts*qcbars/max(D_dts,Dmin) else if ( D_dts.gt.Dmin ) then tmp2s = D_dts*qcbars/max(D_dts,Dmin) else tmp2s = 0. endif endif !do phase changes from large-scale processes ST(i,k,j) = ST(i,k,j) + hls*max(dcond_ls_ice(i,k),0.)/cp_air - & hls*(max(-1.*dcond_ls_ice(i,k),0.)/max(qi_mean(i,k),qmin))*tmp2s/cp_air SQ(i,k,j) = SQ(i,k,j) - max(dcond_ls_ice(i,k),0.) + & (max(-1.*dcond_ls_ice(i,k),0.)/max(qi_mean(i,k),qmin))*tmp2s !cloud erosion changes temperature and vapor ST(i,k,j) = ST(i,k,j) - hls*D_eros(i,k)* tmp2s/cp_air SQ(i,k,j) = SQ(i,k,j) + D_eros(i,k)* tmp2s !add settling ice flux to snow_cld snow_cld(i,k) = D1_dt(i,k)*tmp2s*deltpg(i,k)*inv_dtcloud !add melting of ice to temperature tendency ST(i,k,j) = ST(i,k,j) - hlf*D2_dt(i,k)*tmp2s/cp_air !add melting of ice to the rainflux rain_cld(i,k) = rain_cld(i,k) + D2_dt(i,k)*tmp2s*deltpg(i,k)*inv_dtcloud ! ! diagnostics for cloud ice tendencies ! ! if (max(diag_id%qidt_dep,diag_id%qi_dep_col) > 0) & ! qidt_dep (i,k,j) = max(dcond_ls_ice(i,k),0.)*inv_dtcloud if (max(diag_id%qidt_dep,diag_id%qi_dep_col) > 0) & diag_4d(i,k,j,diag_pt%qidt_dep ) = max(dcond_ls_ice(i,k),0.)*inv_dtcloud ! if (max(diag_id%qidt_subl,diag_id%qi_subl_col) > 0) & ! qidt_subl(i,k,j) = (max(0.,-1.*dcond_ls_ice(i,k))/ & ! max(qi_mean(i,k),qmin))*tmp2s*inv_dtcloud if (max(diag_id%qidt_subl,diag_id%qi_subl_col) > 0) & diag_4d(i,k,j,diag_pt%qidt_subl) = (max(0.,-1.*dcond_ls_ice(i,k))/ & max(qi_mean(i,k),qmin))*tmp2s*inv_dtcloud ! if (max(diag_id%qidt_melt,diag_id%qi_melt_col) > 0) & ! qidt_melt(i,k,j) = D2_dt(i,k) *tmp2s*inv_dtcloud if (max(diag_id%qidt_melt,diag_id%qi_melt_col) > 0) & diag_4d(i,k,j,diag_pt%qidt_melt) = D2_dt(i,k) *tmp2s*inv_dtcloud ! if (max(diag_id%qidt_eros,diag_id%qi_eros_col) > 0) & ! qidt_eros(i,k,j) = D_eros(i,k)*tmp2s*inv_dtcloud if (max(diag_id%qidt_eros,diag_id%qi_eros_col) > 0) & diag_4d(i,k,j,diag_pt%qidt_eros) = D_eros(i,k)*tmp2s*inv_dtcloud if (max(diag_id%lsf_strat,diag_id%lcf_strat,diag_id%mfls_strat) > 0) then tmp1s = max(dcond_ls_ice(i,k),0.)*inv_dtcloud !RSH BUGFIX 6/25/10: ! if ( (qldt_cond(i,k,j) + tmp1s) .gt. 0. ) lsf_strat(:,:,j) = 1. ! if ( (qldt_cond(i,k,j) + tmp1s) .gt. 0. ) diag_4d(:,:,j, diag_pt%lsf_strat) = 1. ! if ( (qldt_cond(i,k,j) + tmp1s) .gt. 0. ) lsf_strat(i,k,j) = 1. ! if ( (qldt_cond(i,k,j) + tmp1s) .gt. 0. ) diag_4d(i,k,j, diag_pt%lsf_strat) = 1. ! if ( (diag_4d(i,k,j,diag_pt%qldt_cond) + tmp1s) .gt. 0. ) lsf_strat(i,k,j) = 1. if ( (diag_4d(i,k,j,diag_pt%qldt_cond) + tmp1s) .gt. 0. ) diag_4d(i,k,j, diag_pt%lsf_strat) = 1. endif enddo enddo !----- -----! ! ! ! END OF ! ! ! ! ICE PHASE MICROPHYSICS ! ! ! ! AND ! ! ! ! ANALYTIC INTEGRATION OF QI EQUATION ! ! ! ! SECTION ! ! ! ! ! ! ! !----------------------------------------------------------------------! !----------------------------------------------------------------------- ! ! ! ! RAIN EVAPORATION ! ! ! Rain evaporation is derived by integration of the growth ! equation of a droplet over the assumed Marshall-Palmer ! distribution of rain drops (equation #22). This leads to the ! following formula: ! ! (50) dqv/dt_local = 56788.636 * {rain_rate/dens_h2o}^(11/18) *(1-U)/ ! ! ( SQRT(airdens)* A_plus_B) ! ! Numerically this equation integrated by use of time-centered ! values for qs and qv. This leads to the solution: ! ! (51) qv_clr(t+1)-qv_clr(t) = K3 *[qs(t)-qv_clr(t)]/ ! {1.+0.5*K3*(1+gamma)} ! where ! ! K3= 56788.636 * dtcloud * {rain_rate_local/dens_h2o}^(11/18) / ! ( SQRT(airdens)* A_plus_B * qs) ! ! and gamma is given by (3). Note that in (51), it is made ! explicit that it is the vapor concentration in the unsaturated ! part of the grid box that is used in the rain evaporation ! formula. ! ! Now there are several limiters to this formula. First, ! you cannot evaporate more than is available in a time step. ! The amount available for evaporation locally is ! (rain_clr/a_rain_clr)*(grav*dtcloud/deltp). Second, to ! avoid supersaturating the box or exceeding the critical ! relative humidity above which rain does not evaporate, ! the amount of evaporation is limited. ! ! Finally rain evaporation occurs only if the relative humidity ! in the unsaturated portion of the grid box, U_clr, is less ! then a threshold, U_evap. U_evap, will not necessarily be ! one. For example, stratiform precipitation in convective ! regions rarely saturates subcloud air because of the presence ! of downdrafts. If the convection scheme does not have down- ! drafts then it doesn't make sense to allow the sub-cloud layer ! to saturate. U_clr may be solved from (8) as: ! ! (52) U_clr = ( U - qa ) / (1. - qa) ! ! Some variables are temporarily stored in tmp1. ! ! Note that for pdf clouds the relative humidity in the clear part ! of the grid box can be calculated exactly from the beta distr- ! ibution. do k=1,jdim do i=1,idim !compute U_clr if (.not. do_pdf_clouds) then U_clr(i,k) = (U(i,k)-qa_mean(i,k))/max((1.-qa_mean(i,k)),qmin) else U_clr(i,k) = qvg(i,k)/qs(i,k,j) end if !keep U_clr > 0. and U_clr < 1. U_clr(i,k) = min(max(U_clr(i,k),0.),1.) !compute K3 tmp1s = 56788.636 * dtcloud * ((rain_clr(i,k)/max(a_rain_clr(i,k),qmin)/ & dens_h2o)**(11./18.))/SQRT(airdens(i,k,j))/A_plus_B(i,k,j)& /qs(i,k,j) !compute local change in vapor mixing ratio due to !rain evaporation tmp1s = tmp1s*qs(i,k,j)*(1.-U_clr(i,k))/(1.+0.5*tmp1s*(1.+gamma(i,k,j))) !limit change in qv to the amount that would raise the relative !humidity to U_evap in the clear portion of the grid box tmp1s = min(tmp1s,((1.-qa_mean(i,k))/max(a_rain_clr(i,k),qmin))*qs(i,k,j)* & max(0.,U_evap-U_clr(i,k))/(1.+(U_evap*(1.-qa_mean(i,k))+qa_mean(i,k))* & gamma(i,k,j)) ) !do limiter by amount available tmp1s= tmp1s*a_rain_clr(i,k)*deltpg(i,k)*inv_dtcloud tmp2s= max(min(rain_clr(i,k),tmp1s),0.) SQ(i,k,j) = SQ(i,k,j) + tmp2s*dtcloud/deltpg(i,k) ST(i,k,j) = ST(i,k,j) - hlv*tmp2s*dtcloud/deltpg(i,k)/cp_air !if all of the rain evaporates set things to zero. if (tmp1s.gt.rain_clr(i,k).and.a_rain_clr(i,k).gt.qmin) then rain_clr(i,k) = 0. a_rain_clr(i,k) = 0. else rain_clr(i,k) = rain_clr(i,k) - tmp2s endif ! if (max(diag_id%rain_evap,diag_id%rain_evap_col) > 0) rain_evap(i,k,j) = tmp2s/deltpg(i,k) if (max(diag_id%rain_evap,diag_id%rain_evap_col) > 0) diag_4d(i,k,j,diag_pt%rain_evap) = tmp2s/deltpg(i,k) !rab enddo !rab enddo !----------------------------------------------------------------------- ! ! ! ! SNOW SUBLIMATION ! ! ! Sublimation of cloud ice ! ! [The following follows Rotstayn (1997)] ! Given the assumptions of the Marshall-Palmer distribution of ! ice crystals (18), the crystal growth equation as a function ! of the humidity of the air and the diffusivity of water vapor ! and thermal conductivity of air is integrated over all crystal ! sizes. This yields: ! ! (53) dqi/dt_local = - a_snow_clr* K3 * (qs - qv_clr) ! ! where the leading factor of a_snow_clr is the portion of the ! grid box undergoing sublimation. K3 is given by ! ! (54) K3 = (4/(pi*rho_air*qs*rho_ice*A_plus_B))* ! ((snow_clr/a_snow_clr/3.29)**1.16 ) * ! [ 0.65*lamda_f^2 + ! 198.92227 * (airdens)^0.5 * ! ((snow_clr/a_snow_clr)**(1/14.5)) * lamda_f^(3/2) ] ! ! Note that equation (53) is identical to equation (30) of ! Rotstayn. ! ! Numerically this is integrated as in rain evaporation. !rab do k=1,jdim !rab do i=1,idim !compute K3 tmp1s = dtcloud * (4./3.14159/rho_ice/airdens(i,k,j)/ & A_plus_B(i,k,j)/qs(i,k,j))*((snow_clr(i,k)/max(a_snow_clr(i,k), & qmin)/3.29)**(1./1.16))*(0.65*lamda_f(i,k)*lamda_f(i,k) + & 198.92227*lamda_f(i,k)*SQRT(airdens(i,k,j)*lamda_f(i,k))* & ( (snow_clr(i,k)/max(a_snow_clr(i,k),qmin))**(1./14.5) ) ) !compute local change in vapor mixing ratio due to !snow sublimation tmp1s = tmp1s*qs(i,k,j)*(1.-U_clr(i,k))/(1.+0.5*tmp1s*(1.+gamma(i,k,j))) !limit change in qv to the amount that would raise the relative !humidity to U_evap in the clear portion of the grid box tmp1s = min(tmp1s,((1.-qa_mean(i,k))/max(a_snow_clr(i,k),qmin))*qs(i,k,j)* & max(0.,U_evap-U_clr(i,k))/(1.+(U_evap*(1.-qa_mean(i,k))+qa_mean(i,k))* & gamma(i,k,j)) ) !do limiter by amount available tmp1s= tmp1s*a_snow_clr(i,k)*deltpg(i,k)*inv_dtcloud tmp2s= max(min(snow_clr(i,k),tmp1s),0.) SQ(i,k,j) = SQ(i,k,j) + tmp2s*dtcloud/deltpg(i,k) ST(i,k,j) = ST(i,k,j) - hls*tmp2s*dtcloud/deltpg(i,k)/cp_air !if all of the snow sublimates set things to zero. if (tmp1s.gt.snow_clr(i,k).and.a_snow_clr(i,k).gt.qmin) then snow_clr(i,k) = 0. a_snow_clr(i,k) = 0. else snow_clr(i,k) = snow_clr(i,k) - tmp2s endif ! if (max(diag_id%snow_subl,diag_id%snow_subl_col) > 0) snow_subl(i,k,j) = tmp2s/deltpg(i,k) if (max(diag_id%qdt_snow_sublim,diag_id%q_snow_sublim_col) > 0) diag_4d(i,k,j,diag_pt%qdt_snow_sublim) = tmp2s/deltpg(i,k) enddo enddo !----------------------------------------------------------------------- ! ! Adjustment to try to prevent negative water vapor. Adjustment ! will not remove more condensate than available. Cloud amount ! is not adjusted. This is left for the remainder of the ! desctruction code. (cjg) ! if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) qadt_destr(:,:,j) = qadt_destr(:,:,j) + SA(:,:,j)*inv_dtcloud if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) diag_4d(:,:,j,diag_pt%qadt_destr) = diag_4d(:,:,j,diag_pt%qadt_destr) + SA(:,:,j)*inv_dtcloud ! if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) qldt_destr(:,:,j) = qldt_destr(:,:,j) + SL(:,:,j)*inv_dtcloud if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) diag_4d(:,:,j,diag_pt%qldt_destr) = diag_4d(:,:,j,diag_pt%qldt_destr) + SL(:,:,j)*inv_dtcloud ! if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) qidt_destr(:,:,j) = qidt_destr(:,:,j) + SI(:,:,j)*inv_dtcloud if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) diag_4d(:,:,j,diag_pt%qidt_destr) = diag_4d(:,:,j,diag_pt%qidt_destr) + SI(:,:,j)*inv_dtcloud ! if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) qndt_destr(:,:,j) = qndt_destr(:,:,j) + SN(:,:,j)*inv_dtcloud if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) diag_4d(:,:,j,diag_pt%qndt_destr) = diag_4d(:,:,j,diag_pt%qndt_destr) + SN(:,:,j)*inv_dtcloud do k=1,jdim do i=1,idim tmp1s = qv(i,k,j) + SQ(i,k,j) tmp2s = 0.0 tmp3s = 0.0 if ( tmp1s.lt.0.0 ) then if (T(i,k,j).le.tfreeze-40.) then tmp2s = min( -tmp1s, ql_upd(i,k) ) ! liquid to evaporate tmp3s = min( -tmp1s-tmp2s, qi_upd(i,k) ) ! ice to sublimate else tmp3s = min( -tmp1s, qi_upd(i,k) ) ! ice to sublimate tmp2s = min( -tmp1s-tmp3s, ql_upd(i,k) ) ! liquid to evaporate end if ql_upd(i,k) = ql_upd(i,k) - tmp2s qi_upd(i,k) = qi_upd(i,k) - tmp3s SL(i,k,j) = SL(i,k,j) - tmp2s SI(i,k,j) = SI(i,k,j) - tmp3s SQ(i,k,j) = SQ(i,k,j) + tmp2s + tmp3s ST(i,k,j) = ST(i,k,j) - hlv*tmp2s/cp_air - hls*tmp3s/cp_air end if enddo enddo !----------------------------------------------------------------------- ! Cloud Destruction occurs where both ql and qi are .le. qmin, ! or if qa is .le. qmin. In this case, ql, qi, and qa are set to ! zero conserving moisture and energy. if (.not.do_liq_num) then do k=1,jdim do i=1,idim if ((ql_upd(i,k) .le. qmin .and. qi_upd(i,k) .le. qmin) & .or. (qa_upd(i,k) .le. qmin)) then SL(i,k,j) = SL(i,k,j) - ql_upd(i,k) SI(i,k,j) = SI(i,k,j) - qi_upd(i,k) SQ(i,k,j) = SQ(i,k,j) + ql_upd(i,k) + qi_upd(i,k) ST(i,k,j) = ST(i,k,j) - (hlv*ql_upd(i,k) + hls*qi_upd(i,k))/cp_air SA(i,k,j) = SA(i,k,j) - qa_upd(i,k) ql_upd(i,k) = 0.0 qi_upd(i,k) = 0.0 qa_upd(i,k) = 0.0 endif enddo enddo else do k=1,jdim do i=1,idim if ((ql_upd(i,k) .le. qmin .and. qi_upd(i,k) .le. qmin) & .or. (qa_upd(i,k) .le. qmin)) then SL(i,k,j) = SL(i,k,j) - ql_upd(i,k) SI(i,k,j) = SI(i,k,j) - qi_upd(i,k) SQ(i,k,j) = SQ(i,k,j) + ql_upd(i,k) + qi_upd(i,k) ST(i,k,j) = ST(i,k,j) - (hlv*ql_upd(i,k) + hls*qi_upd(i,k))/cp_air SA(i,k,j) = SA(i,k,j) - qa_upd(i,k) SN(i,k,j) = SN(i,k,j) - qn_upd(i,k) ql_upd(i,k) = 0.0 qi_upd(i,k) = 0.0 qa_upd(i,k) = 0.0 qn_upd(i,k) = 0.0 endif enddo enddo endif ! if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) qadt_destr(:,:,j) = qadt_destr(:,:,j) - SA(:,:,j)*inv_dtcloud if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) diag_4d(:,:,j,diag_pt%qadt_destr) = diag_4d(:,:,j,diag_pt%qadt_destr) - SA(:,:,j)*inv_dtcloud ! if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) qldt_destr(:,:,j) = qldt_destr(:,:,j) - SL(:,:,j)*inv_dtcloud if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) diag_4d(:,:,j,diag_pt%qldt_destr) = diag_4d(:,:,j,diag_pt%qldt_destr) - SL(:,:,j)*inv_dtcloud ! if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) qidt_destr(:,:,j) = qidt_destr(:,:,j) - SI(:,:,j)*inv_dtcloud if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) diag_4d(:,:,j,diag_pt%qidt_destr) = diag_4d(:,:,j,diag_pt%qidt_destr) - SI(:,:,j)*inv_dtcloud ! if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) qndt_destr(:,:,j) = qndt_destr(:,:,j) - SN(:,:,j)*inv_dtcloud if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) diag_4d(:,:,j,diag_pt%qndt_destr) = diag_4d(:,:,j,diag_pt%qndt_destr) - SN(:,:,j)*inv_dtcloud !----------------------------------------------------------------------- ! ! Adjustment. Due to numerical errors in detrainment or advection ! sometimes the current state of the grid box may be super- ! saturated. Under the assumption that the temperature is constant ! in the grid box and that q <= qs, the excess vapor is condensed. ! ! What happens to the condensed water vapor is determined by the ! namelist parameter super_choice. ! ! If super_choice = .false. (default), then the condensed water is ! is added to the precipitation fluxes. If super_choice = .true., ! then the condensed water is added to the cloud condensate field. ! Note that in this case the cloud fraction is raised to one. ! ! The phase partitioning depends on super_choice; if super_choice ! is false then at T < -20C, snow is produced. If super_choice ! is true, then at T < -40C, ice is produced. The latter choice ! is consistent with that done in the large-scale condensation ! section above. ! ! If pdf clouds are operating then this section is bypassed - ! as statistical clouds should remove supersaturation according ! to the beta distribution used. if (.not.do_pdf_clouds) then ! Old code ! !estimate current qs ! tmp2 = qs(:,:,j)+dqsdT(:,:,j)*ST(:,:,j) ! ! !compute excess over saturation ! tmp1 = max(0.,qv(:,:,j)+SQ(:,:,j)-tmp2)/(1.+gamma(:,:,j)) ! New more accurate version !RSH 9/18/09: should be within the super_choice if block: ! updated temperature in tmp1 tmp1 = T(:,:,j) + ST(:,:,j) ! updated qs in tmp2, updated gamma in tmp3, updated dqsdT in tmp5 call compute_qs(tmp1, pfull(:,:,j), tmp2, dqsdT=tmp5) tmp3 = tmp5 *(min(1.,max(0.,0.05*(tmp1-tfreeze+20.)))*hlv + & min(1.,max(0.,0.05*(tfreeze -tmp1 )))*hls)/cp_air !compute excess over saturation tmp1 = max(0.,qv(:,:,j)+SQ(:,:,j)-tmp2)/(1.+tmp3) !rab - save off tmp1 for diagnostic ice/liq adjustment fields in liq_adj array ! if (max(diag_id%ice_adj,diag_id%ice_adj_col,diag_id%liq_adj,diag_id%liq_adj_col) .gt. 0) & ! liq_adj(:,:,j) = tmp1*inv_dtcloud if (max(diag_id%ice_adj,diag_id%ice_adj_col,diag_id%liq_adj,diag_id%liq_adj_col) .gt. 0) & diag_4d(:,:,j,diag_pt%liq_adj) = tmp1*inv_dtcloud !change vapor content SQ(:,:,j)=SQ(:,:,j)-tmp1 if (super_choice) then ! Put supersaturation into cloud !cloud fraction source diagnostic if (max(diag_id%qadt_super,diag_id%qa_super_col) > 0) then where (tmp1 .gt. 0.) ! qadt_super(:,:,j) = (1.-qa_upd) * inv_dtcloud diag_4d(:,:,j,diag_pt%qadt_super) = (1.-qa_upd) * inv_dtcloud endwhere endif !yim 11/7/07 if (do_liq_num) then do k=1,jdim do i=1,idim if (T(i,k,j) > tfreeze - 40. .and. & tmp1(i,k) > 0.0) then qn_upd(i,k) = qn_upd(i,k) + drop1(i,k)*1.0e6/ & airdens(i,k,j)*(1. - qa_upd(i,k)) SN(i,k,j) = SN(i,k,j) + drop1(i,k)*1.e6/ & airdens(i,k,j)*(1.-qa_upd(i,k)) if (max(diag_id%qndt_super,diag_id%qn_super_col) > 0) & ! qndt_super(i,k,j) = qndt_super(i,k,j) + & diag_4d(i,k,j,diag_pt%qndt_super) = diag_4d(i,k,j,diag_pt%qndt_super) + & drop1(i,k)*1.e6 / & airdens(i,k,j)*(1.-qa_upd(i,k))*inv_dtcloud endif end do end do ! if (max(diag_id%qndt_super,diag_id%qn_super_col) > 0) & ! diag_4d(:,:,j,diag_pt%qndt_super) = qndt_super(:,:,j) endif !add in excess to cloud condensate, change cloud area and !increment temperature do k=1,jdim do i=1,idim if(tmp1(i,k).gt.0) then if (T(i,k,j) .le. tfreeze-40.)then qi_upd(i,k) = qi_upd(i,k) + tmp1(i,k) SI(i,k,j) = SI(i,k,j) + tmp1(i,k) ST(i,k,j) = ST(i,k,j) + hls*tmp1(i,k)/cp_air else ! where (T(i,k,j) .gt. tfreeze-40.) ql_upd(i,k) = ql_upd(i,k) + tmp1(i,k) SL(i,k,j) = SL(i,k,j) + tmp1(i,k) ST(i,k,j) = ST(i,k,j) + hlv*tmp1(i,k)/cp_air endif if (limit_conv_cloud_frac) then tmp2s = ahuco(i,k,j) else tmp2s = 0. endif SA(i,k,j) = SA(i,k,j) + (1.-qa_upd(i,k)-tmp2s) qa_upd(i,k) = 1. - tmp2s endif enddo enddo else !Put supersaturation into precip !add in excess to precipitation fluxes, change their area !and increment temperature do k=1,jdim do i=1,idim if(tmp1(i,k).gt.0) then if (T(i,k,j) .le. tfreeze-20.) then snow_cld(i,k) = snow_cld(i,k) + qa_mean(i,k) *tmp1(i,k)*deltpg(i,k)*inv_dtcloud snow_clr(i,k) = snow_clr(i,k) + (1.-qa_mean(i,k))*tmp1(i,k)*deltpg(i,k)* & inv_dtcloud a_snow_cld(i,k) = qa_mean(i,k) a_snow_clr(i,k) = 1.-qa_mean(i,k) ST(i,k,j) = ST(i,k,j) + hls*tmp1(i,k)/cp_air else ! where (T(i,k,j) .gt. tfreeze-20.) rain_cld(i,k) = rain_cld(i,k) + qa_mean(i,k) *tmp1(i,k)*deltpg(i,k)*inv_dtcloud rain_clr(i,k) = rain_clr(i,k) + (1.-qa_mean(i,k))*tmp1(i,k)*deltpg(i,k)* & inv_dtcloud a_rain_cld(i,k) = qa_mean(i,k) a_rain_clr(i,k) = 1.-qa_mean(i,k) ST(i,k,j) = ST(i,k,j) + hlv*tmp1(i,k)/cp_air endif endif enddo enddo end if !super choice end if !for do_pdf_clouds !----------------------------------------------------------------------- ! Final clean up to remove numerical noise ! ! if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) qadt_destr(:,:,j) = qadt_destr(:,:,j) + SA(:,:,j)*inv_dtcloud if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) diag_4d(:,:,j,diag_pt%qadt_destr) = diag_4d(:,:,j,diag_pt%qadt_destr) + SA(:,:,j)*inv_dtcloud ! if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) qldt_destr(:,:,j) = qldt_destr(:,:,j) + SL(:,:,j)*inv_dtcloud if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) diag_4d(:,:,j,diag_pt%qldt_destr) = diag_4d(:,:,j,diag_pt%qldt_destr) + SL(:,:,j)*inv_dtcloud ! if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) qidt_destr(:,:,j) = qidt_destr(:,:,j) + SI(:,:,j)*inv_dtcloud if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) diag_4d(:,:,j,diag_pt%qidt_destr) = diag_4d(:,:,j,diag_pt%qidt_destr) + SI(:,:,j)*inv_dtcloud ! if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) qndt_destr(:,:,j) = qndt_destr(:,:,j) + SN(:,:,j)*inv_dtcloud if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) diag_4d(:,:,j,diag_pt%qndt_destr) = diag_4d(:,:,j,diag_pt%qndt_destr) + SN(:,:,j)*inv_dtcloud do k=1,jdim do i=1,idim ql_upd(i,k) = ql(i,k,j) + SL(i,k,j) if ( abs(ql_upd(i,k)) .le. qmin & .and. qv(i,k,j)+SQ(i,k,j)+ql_upd(i,k) > 0.0 ) then SL(i,k,j) = -ql(i,k,j) SQ(i,k,j) = SQ(i,k,j) + ql_upd(i,k) ST(i,k,j) = ST(i,k,j) - hlv*ql_upd(i,k)/cp_air ql_upd(i,k) = 0.0 endif qi_upd(i,k) = qi(i,k,j) + SI(i,k,j) if ( abs(qi_upd(i,k)) .le. qmin & .and. qv(i,k,j)+SQ(i,k,j)+qi_upd(i,k) > 0.0 ) then SI(i,k,j) = -qi(i,k,j) SQ(i,k,j) = SQ(i,k,j) + qi_upd(i,k) ST(i,k,j) = ST(i,k,j) - hls*qi_upd(i,k)/cp_air qi_upd(i,k) = 0.0 endif qa_upd(i,k) = qa(i,k,j) + SA(i,k,j) if ( abs(qa_upd(i,k)) .le. qmin ) then SA(i,k,j) = -qa(i,k,j) qa_upd(i,k) = 0.0 endif enddo enddo ! if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) qadt_destr(:,:,j) = qadt_destr(:,:,j) - SA(:,:,j)*inv_dtcloud if (max(diag_id%qadt_destr,diag_id%qa_destr_col) > 0) diag_4d(:,:,j,diag_pt%qadt_destr) = diag_4d(:,:,j,diag_pt%qadt_destr) - SA(:,:,j)*inv_dtcloud ! if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) qldt_destr(:,:,j) = qldt_destr(:,:,j) - SL(:,:,j)*inv_dtcloud if (max(diag_id%qldt_destr,diag_id%ql_destr_col) > 0) diag_4d(:,:,j,diag_pt%qldt_destr) = diag_4d(:,:,j,diag_pt%qldt_destr) - SL(:,:,j)*inv_dtcloud ! if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) qidt_destr(:,:,j) = qidt_destr(:,:,j) - SI(:,:,j)*inv_dtcloud if (max(diag_id%qidt_destr,diag_id%qi_destr_col) > 0) diag_4d(:,:,j,diag_pt%qidt_destr) = diag_4d(:,:,j,diag_pt%qidt_destr) - SI(:,:,j)*inv_dtcloud ! if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) qndt_destr(:,:,j) = qndt_destr(:,:,j) - SN(:,:,j)*inv_dtcloud if (max(diag_id%qndt_destr,diag_id%qn_destr_col) > 0) diag_4d(:,:,j,diag_pt%qndt_destr) = diag_4d(:,:,j,diag_pt%qndt_destr) - SN(:,:,j)*inv_dtcloud !----------------------------------------------------------------------- ! ! Save qa_mean of current level into qa_mean_lst. This is used ! in transferring rain and snow fluxes between levels. qa_mean_lst = qa_mean !----------------------------------------------------------------------- ! ! add the ice falling out from cloud to qidt_fall ! if (max(diag_id%qidt_fall,diag_id%qi_fall_col) > 0) qidt_fall(:,:,j) = qidt_fall(:,:,j) + & ! (snow_cld/deltpg) if (max(diag_id%qidt_fall,diag_id%qi_fall_col) > 0) diag_4d(:,:,j,diag_pt%qidt_fall) = diag_4d(:,:,j,diag_pt%qidt_fall) + & (snow_cld/deltpg) !----------------------------------------------------------------------- ! ! save profiles of rain and snow ! rain3d(:,:,j+1) = rain_clr(:,:) + rain_cld(:,:) snow3d(:,:,j+1) = snow_clr(:,:) + snow_cld(:,:) snowclr3d(:,:,j+1) = snow_clr(:,:) !----------------------------------------------------------------------- ! ! Save rain and snow diagnostics ! if (diag_id%rain_clr > 0) rain_clr_diag(:,:,j+1) = rain_clr ! if (diag_id%rain_cld > 0) rain_cld_diag(:,:,j+1) = rain_cld ! if (diag_id%a_rain_clr > 0) a_rain_clr_diag(:,:,j+1) = a_rain_clr ! if (diag_id%a_rain_cld > 0) a_rain_cld_diag(:,:,j+1) = a_rain_cld ! if (diag_id%snow_clr > 0) snow_clr_diag(:,:,j+1) = snow_clr ! if (diag_id%snow_cld > 0) snow_cld_diag(:,:,j+1) = snow_cld ! if (diag_id%a_snow_clr > 0) a_snow_clr_diag(:,:,j+1) = a_snow_clr ! if (diag_id%a_snow_cld > 0) a_snow_cld_diag(:,:,j+1) = a_snow_cld ! if (diag_id%a_precip_clr > 0) a_precip_clr_diag(:,:,j+1) = max(a_rain_clr,a_snow_clr) ! if (diag_id%a_precip_cld > 0) a_precip_cld_diag(:,:,j+1) = max(a_rain_cld,a_snow_cld) if (diag_id%rain_clr > 0) diag_4d_kp1(:,:,j+1,diag_pt%rain_clr) = rain_clr if (diag_id%rain_cld > 0) diag_4d_kp1(:,:,j+1,diag_pt%rain_cld) = rain_cld if (diag_id%a_rain_clr > 0) diag_4d_kp1(:,:,j+1,diag_pt%a_rain_clr) = a_rain_clr if (diag_id%a_rain_cld > 0) diag_4d_kp1(:,:,j+1,diag_pt%a_rain_cld) = a_rain_cld if (diag_id%snow_clr > 0) diag_4d_kp1(:,:,j+1,diag_pt%snow_clr) = snow_clr if (diag_id%snow_cld > 0) diag_4d_kp1(:,:,j+1,diag_pt%snow_cld) = snow_cld if (diag_id%a_snow_clr > 0) diag_4d_kp1(:,:,j+1,diag_pt%a_snow_clr) = a_snow_clr if (diag_id%a_snow_cld > 0) diag_4d_kp1(:,:,j+1,diag_pt%a_snow_cld) = a_snow_cld if (diag_id%a_precip_clr > 0) diag_4d_kp1(:,:,j+1,diag_pt%a_precip_clr) = max(a_rain_clr,a_snow_clr) if (diag_id%a_precip_cld > 0) diag_4d_kp1(:,:,j+1,diag_pt%a_precip_cld) = max(a_rain_cld,a_snow_cld) !----------------------------------------------------------------------- ! ! ! Put rain and ice fluxes into surfrain and surfsnow if the ! grid point is at the bottom of a column. If MASK is not ! present then this code is executed only if j .eq. kdim. ! IF MASK is present some grid points may be beneath ground. ! If a given grid point is at the bottom of the column then ! the surface values of rain and snow must be created. ! Also if the MASK is present then the code forces all tenden- ! cies below ground to be zero. Note that MASK = 1. equals above ! ground point, MASK = 0. equals below ground point. if (present(MASK)) then !zero out all tendencies below ground ST(:,:,j)=MASK(:,:,j)*ST(:,:,j) SQ(:,:,j)=MASK(:,:,j)*SQ(:,:,j) SL(:,:,j)=MASK(:,:,j)*SL(:,:,j) SI(:,:,j)=MASK(:,:,j)*SI(:,:,j) SA(:,:,j)=MASK(:,:,j)*SA(:,:,j) if (do_liq_num) then SN(:,:,j)=MASK(:,:,j)*SN(:,:,j) endif if (j .lt. kdim) then !bottom of true points in columns which contain some !dummy points where(MASK(:,:,j) .eq. 1. .and. MASK(:,:,j+1) .eq. 0.) surfrain = dtcloud*(rain_clr+rain_cld) surfsnow = dtcloud*(snow_clr+snow_cld) rain_clr = 0. rain_cld = 0. snow_clr = 0. snow_cld = 0. a_rain_clr = 0. a_rain_cld = 0. a_snow_clr = 0. a_snow_cld = 0. end where else !bottom of column for those columns which contain no !dummy points where(MASK(:,:,j) .eq. 1.) surfrain = dtcloud*(rain_clr+rain_cld) surfsnow = dtcloud*(snow_clr+snow_cld) rain_clr = 0. rain_cld = 0. snow_clr = 0. snow_cld = 0. a_rain_clr = 0. a_rain_cld = 0. a_snow_clr = 0. a_snow_cld = 0. end where end if else !do code if we are at bottom of column if (j .eq. kdim) then surfrain = dtcloud*(rain_clr+rain_cld) surfsnow = dtcloud*(snow_clr+snow_cld) end if end if !----------------------------------------------------------------------- ! ! END LOOP OVER VERTICAL LEVELS ! enddo call mpp_clock_end(sc_loop) call mpp_clock_begin(sc_post_loop) !----------------------------------------------------------------------- ! ! INSTANTANEOUS OUTPUT DIAGNOSTICS ! if (num_strat_pts > 0) then do nn=1,num_strat_pts if (strat_pts(1,nn) >= is .and. strat_pts(1,nn) <= ie .and. & strat_pts(2,nn) >= js .and. strat_pts(2,nn) <= je) then ipt=strat_pts(1,nn); jpt=strat_pts(2,nn) i=ipt-is+1; j=jpt-js+1 unit = open_ieee32_file ('strat.data', action='append') write (unit) ipt,jpt, ql(i,j,:)+SL(i,j,:) write (unit) ipt,jpt, qi(i,j,:)+SI(i,j,:) write (unit) ipt,jpt, qa(i,j,:)+SA(i,j,:) write (unit) ipt,jpt, T(i,j,:)+ST(i,j,:) write (unit) ipt,jpt, qv(i,j,:)+SQ(i,j,:) write (unit) ipt,jpt, pfull(i,j,:) call close_file(unit) endif enddo endif !----------------------------------------------------------------------- ! ! DIAGNOSTICS ! !rab - perform the assignments for ice/liq adjustments diagnostics if (max(diag_id%ice_adj,diag_id%ice_adj_col,diag_id%liq_adj,diag_id%liq_adj_col) .gt. 0) then freeze_pt=40. if (.not. super_choice) freeze_pt=20. where (T .le. tfreeze-freeze_pt) ! ice_adj = liq_adj diag_4d(:,:,:,diag_pt%ice_adj) = diag_4d(:,:,:,diag_pt%liq_adj) ! liq_adj = 0. diag_4d(:,:,:,diag_pt%liq_adj) = 0. endwhere endif if ( diag_id%droplets_wtd > 0 ) then !a diag_4d(:,:,:, diag_pt%droplets_wtd) = N3D*ql diag_4d(:,:,:, diag_pt%droplets_wtd) = & diag_4d(:,:,:,diag_pt%droplets)*ql(:,:,:) endif if ( diag_id%subgrid_w_variance > 0 ) then ! debug2=debug2**0.5 ! diag_4d(:,:,:,diag_pt%debug2_3d) = debug2 diag_4d(:,:,:,diag_pt%subgrid_w_variance) = diag_4d(:,:,:,diag_pt%subgrid_w_variance)**0.5 endif !RSH BUGFIX 6/25/10 ! j is undefined here if ( max(diag_id%lcf_strat,diag_id%mfls_strat) > 0 ) then ! where (omega(:,:,j)+grav*Mc(:,:,j) .lt. 0 .and. lsf_strat(:,:,j) .eq.1) ! lcf_strat(:,:,j) = 1. ! diag_4d(:,:,j,diag_pt%lcf_strat) = lcf_strat(:,:,j) ! mfls_strat(:,:,j) = omega(:,:,j)+grav*Mc(:,:,j) ! diag_4d(:,:,j,diag_pt%mfls_strat) = mfls_strat(:,:,j) ! where (omega(:,:,:)+grav*Mc(:,:,:) .lt. 0 .and. lsf_strat(:,:,:) .eq.1) where (omega(:,:,:)+grav*Mc(:,:,:) .lt. 0 .and. diag_4d(:,:,:,diag_pt%lsf_strat) .eq.1) ! lcf_strat(:,:,:) = 1. ! diag_4d(:,:,:,diag_pt%lcf_strat) = lcf_strat(:,:,:) diag_4d(:,:,:,diag_pt%lcf_strat) = 1.0 ! mfls_strat(:,:,:) = omega(:,:,:)+grav*Mc(:,:,:) ! diag_4d(:,:,:,diag_pt%mfls_strat) = mfls_strat(:,:,:) diag_4d(:,:,:,diag_pt%mfls_strat) = omega(:,:,:)+grav*Mc(:,:,:) end where end if do k = 1, kdim do j=1,jdim do i=1,idim deltpg_3d(i,j,k) = (phalf(i,j,k+1)-phalf(i,j,k))/grav if (present(MASK)) then deltpg_3d(i,j,k)=deltpg_3d(i,j,k)*MASK(i,j,k) endif end do end do end do !-------write out column integrated diagnostics------! !GENERIC FORM: do n8=1, n_diag_4d do j = kdim, 1, -1 ! diag_3d(:,:,diag_pt%qldt_evap) = diag_3d(:,:,diag_pt%qldt_evap) & ! + diag_4d(:,:,j,diag_pt%qldt_evap)*deltpg_3d(:,:,j) diag_3d(:,:,n8) = diag_3d(:,:,n8) & + diag_4d(:,:,j,n8)*deltpg_3d(:,:,j) enddo enddo !SPECIAL CASES: !yim: in-cloud droplet column burden if (diag_id%droplets_col > 0) then if (present (qn)) then if (do_liq_num ) then N3D_col(:,:) = 0. do k = 1, kdim do j=1,jdim do i=1,idim ! deltpg(i,j) = (phalf(i,j,k+1)-phalf(i,j,k))/grav ! if (present(MASK)) then ! deltpg(i,j)=deltpg(i,j)*MASK(i,j,k) ! endif if (ql(i,j,k) > qmin .and. & qa(i,j,k) > qmin .and. & qn(i,j,k) > qmin ) then if (qa(i,j,k) > 0.05) then N3D_col(i,j) = N3D_col(i,j) + qn(i,j,k)* & ! airdens(i,j,k)*deltpg(i,j)* & !RSH 12/22/11 fix as per email from yim 11/3/11: ! airdens(i,j,k)*deltpg_3d(i,j,k)* & ! 1.e-6/min(qa(i,j,k),1.) !RSH 12/22/11 !NOTE still differs from new strat_cloud code in that is appplied to ! input fields rather than output fields deltpg_3d(i,j,k)* & 1.e-4/min(qa(i,j,k),1.) endif endif end do end do end do diag_3d(:,:,diag_pt%droplets_col) = N3D_col endif endif endif call mpp_clock_end(sc_post_loop) !----------------------------------------------------------------------- ! ! ! end of subroutine end subroutine strat_cloud_legacy !##################################################################### subroutine strat_cloud_legacy_end (do_pdf_clouds) logical, intent(in) :: do_pdf_clouds if (do_pdf_clouds) call beta_dist_end module_is_initialized = .false. end subroutine strat_cloud_legacy_end !##################################################################### subroutine aerosol_effects (is, js, Time, phalf, airdens, T, diag_4d, & diag_id, diag_pt, & concen_dust_sub, totalmass1, Aerosol, mask) type(diag_id_type), intent(in) :: diag_id type(diag_pt_type), intent(in) :: diag_pt integer, intent (in) :: is,js type(time_type), intent (in) :: Time real, dimension(:,:,:), intent(in ) :: phalf, airdens, T real, dimension(:,:,:,0:), intent(inout ) :: diag_4d real, dimension(:,:,:), intent(out) :: concen_dust_sub real, dimension(:,:,:,:), intent(out) :: totalmass1 type(aerosol_type), intent (in), optional :: Aerosol real, intent (in), optional, dimension(:,:,:) :: mask real, dimension(size(T,1),size(T,2),size(T,3)) :: pthickness real, dimension(size(T,1),size(T,2),size(T,3)) :: concen, & concen_all_sub, & concen_ss_sub, concen_ss_sup,& concen_om, concen_na, concen_an integer :: i,j,k, na , s integer :: idim, jdim, kdim logical :: used idim = size(T,1) jdim = size(T,2) kdim = size(T,3) concen_dust_sub(:,:,:) = 0. totalmass1(:,:,:,:) = 0. if (diag_id%sulfate > 0) then concen_an(:,:,:) = 0. concen_na(:,:,:) = 0. concen(:,:,:) = 0. endif if (use_online_aerosol) then concen_ss_sub(:,:,:) = 0. concen_ss_sup(:,:,:) = 0. concen_all_sub(:,:,:) = 0. endif do k = 1,kdim do j = 1,jdim do i = 1,idim if (phalf(i,j,k) < 1.0) then pthickness(i,j,k) = (phalf(i,j,k+1) - phalf(i,j,k))/& grav/airdens(i,j,k) else pthickness(i,j,k) = log(phalf(i,j,k+1)/ & phalf(i,j,k))*8.314*T(i,j,k)/(9.8*0.02888) end if end do end do end do if (present (Aerosol)) then if (do_liq_num) then if (use_online_aerosol) then do na = 1,size(Aerosol%aerosol,4) if (trim(Aerosol%aerosol_names(na)) == 'so4' .or. & trim(Aerosol%aerosol_names(na)) == 'so4_anthro' .or.& trim(Aerosol%aerosol_names(na)) == 'so4_natural') & then do k = 1,kdim do j = 1,jdim do i = 1,idim totalmass1(i,j,k,1) = totalmass1(i,j,k,1) + & Aerosol%aerosol(i,j,k,na) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'omphilic' .or.& trim(Aerosol%aerosol_names(na)) == 'omphobic') & then do k = 1,kdim do j = 1,jdim do i = 1,idim totalmass1(i,j,k,4) = totalmass1(i,j,k,4) + & Aerosol%aerosol(i,j,k,na) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'seasalt1' .or.& trim(Aerosol%aerosol_names(na)) == 'seasalt2') & then do k = 1,kdim do j = 1,jdim do i = 1,idim concen_ss_sub(i,j,k) = concen_ss_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'seasalt3' .or.& trim(Aerosol%aerosol_names(na)) == 'seasalt4' .or.& trim(Aerosol%aerosol_names(na)) == 'seasalt5') & then do k = 1,kdim do j = 1,jdim do i = 1,idim concen_ss_sup(i,j,k) = concen_ss_sup(i,j,k) + & Aerosol%aerosol(i,j,k,na) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'bcphilic' .or.& trim(Aerosol%aerosol_names(na)) == 'bcphobic' .or.& trim(Aerosol%aerosol_names(na)) == 'dust1' .or.& trim(Aerosol%aerosol_names(na)) == 'dust2' .or.& trim(Aerosol%aerosol_names(na)) == 'dust3') & then do k = 1,kdim do j = 1,jdim do i = 1,idim concen_all_sub(i,j,k) = concen_all_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) end do end do end do endif if (do_dust_berg) then if (trim(Aerosol%aerosol_names(na)) == 'dust1' .or. & trim(Aerosol%aerosol_names(na)) == 'dust2' .or. & trim( Aerosol%aerosol_names(na)) == 'dust3') then do k = 1,kdim do j = 1,jdim do i = 1,idim concen_dust_sub(i,j,k) = & concen_dust_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) end do end do end do endif endif end do do k = 1,kdim do j = 1,jdim do i = 1,idim totalmass1(i,j,k,3) = concen_ss_sub(i,j,k) totalmass1(i,j,k,2) = concen_all_sub(i,j,k) + & totalmass1(i,j,k,4) + & concen_ss_sub(i,j,k) end do end do end do if (use_sub_seasalt) then else do k = 1,kdim do j = 1,jdim do i = 1,idim totalmass1(i,j,k,3) = concen_ss_sub(i,j,k) + & concen_ss_sup(i,j,k) end do end do end do endif if (diag_id%sulfate > 0) then do k = 1,kdim do j = 1,jdim do i = 1,idim concen(i,j,k) = 0.7273*totalmass1(i,j,k,1)/ & pthickness(i,j,k)*1.0e9 end do end do end do endif do k = 1,kdim do j = 1,jdim do i = 1,idim concen_ss_sub(i,j,k) = concen_ss_sub(i,j,k)/ & pthickness(i,j,k)*1.0e9 concen_ss_sup(i,j,k) = concen_ss_sup(i,j,k)/ & pthickness(i,j,k)*1.0e9 end do end do end do else ! (use_online_aerosol) if (do_dust_berg) then ! YMice submicron dust (NO. 14 to NO. 18) do s = 14,18 do k = 1,kdim do j = 1,jdim do i = 1,idim concen_dust_sub(i,j,k) = concen_dust_sub(i,j,k)+ & Aerosol%aerosol(i,j,k,s) end do end do end do end do endif if (diag_id%sulfate > 0) then do k = 1,kdim do j = 1,jdim do i = 1,idim ! anthro. and natural sulfate concentration (ug so4/m3) concen_an(i,j,k) = 0.7273*Aerosol%aerosol(i,j,k,1)/& pthickness(i,j,k)*1.0e9 concen_na(i,j,k) = 0.7273*Aerosol%aerosol(i,j,k,2)/& pthickness(i,j,k)*1.0e9 concen(i,j,k) = concen_an(i,j,k) + concen_na(i,j,k) end do end do end do endif do k = 1,kdim do j = 1,jdim do i = 1,idim !offline ! NO. 1 natural Sulfate; NO. 2 anthro. sulfate; NO. 3 Sea Salt; NO. 4 Or ganics totalmass1(i,j,k,1) = Aerosol%aerosol(i,j,k,2) totalmass1(i,j,k,2) = Aerosol%aerosol(i,j,k,1) totalmass1(i,j,k,3) = sea_salt_scale* & Aerosol%aerosol(i,j,k,5) totalmass1(i,j,k,4) = om_to_oc* & Aerosol%aerosol(i,j,k,3) end do end do end do endif ! (use_online_aerosol) do na = 1, 4 do k = 1,kdim do j = 1,jdim do i = 1,idim totalmass1(i,j,k,na) = totalmass1(i,j,k,na)/ & pthickness(i,j,k)*1.0e9*1.0e-12 end do end do end do end do if (do_dust_berg) then do k = 1,kdim do j = 1,jdim do i = 1,idim ! submicron dust concentration (ug/m3) (NO. 2 to NO. 4) concen_dust_sub(i,j,k) = concen_dust_sub(i,j,k)/ & pthickness(i,j,k)*1.0e9 end do end do end do endif if (diag_id%sulfate > 0) then diag_4d(:,:,:,diag_pt%sulfate) = concen (:,:,:) endif if (use_online_aerosol) then if (diag_id%seasalt_sub > 0) then diag_4d(:,:,:,diag_pt%seasalt_sub) = concen_ss_sub endif if (diag_id%seasalt_sup > 0) then diag_4d(:,:,:,diag_pt%seasalt_sup) = concen_ss_sup endif endif if (diag_id%om > 0) then do k = 1,kdim do j = 1,jdim do i = 1,idim concen_om(i,j,k) = totalmass1(i,j,k,2)*1.0e12 diag_4d(i,j,k,diag_pt%om) = concen_om(i,j,k) end do end do end do ! used = send_data (diag_id%om, diag_4d(:,:,:,diag_pt%om), Time, is, js, 1,& ! rmask=mask ) endif endif ! (do_liq_num) endif ! (Present(Aerosol)) !---------------------------------------------------------------------- end subroutine aerosol_effects !####################################################################### !####################################################################### !----------------------------------------------------------------------- !BOP ! !ROUTINE: ppm2m_sak --- Piecewise parabolic method for fields ! ! !INTERFACE: subroutine ppm2m_sak(a4, delp, km, kmap, i1, i2, iv, kord) implicit none ! !INPUT PARAMETERS: integer, intent(in):: iv ! iv =-1: winds ! iv = 0: positive definite scalars ! iv = 1: others integer, intent(in):: i1 ! Starting longitude integer, intent(in):: i2 ! Finishing longitude integer, intent(in):: km ! vertical dimension integer, intent(in):: kmap ! partial remap to start integer, intent(in):: kord ! Order (or more accurately method no.): ! real, intent(in):: delp(i1:i2,km) ! layer pressure thickness ! !INPUT/OUTPUT PARAMETERS: real, intent(inout):: a4(4,i1:i2,km) ! Interpolated values ! !DESCRIPTION: ! ! Perform the piecewise parabolic method ! ! !REVISION HISTORY: ! ??.??.?? Lin Creation ! 02.04.04 Sawyer Newest release from FVGCM ! !EOP !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: ! local arrays: real dc(i1:i2,km) real h2(i1:i2,km) real delq(i1:i2,km) real df2(i1:i2,km) real d4(i1:i2,km) ! local scalars: integer i, k, km1, lmt integer it real fac real a1, a2, c1, c2, c3, d1, d2 real qmax, qmin, cmax, cmin real qm, dq, tmp real qmp, pmp real lac km1 = km - 1 it = i2 - i1 + 1 do k=max(2,kmap-2),km do i=i1,i2 delq(i,k-1) = a4(1,i,k) - a4(1,i,k-1) d4(i,k ) = delp(i,k-1) + delp(i,k) enddo enddo do k=max(2,kmap-2),km1 do i=i1,i2 c1 = (delp(i,k-1)+0.5*delp(i,k))/d4(i,k+1) c2 = (delp(i,k+1)+0.5*delp(i,k))/d4(i,k) tmp = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) / & (d4(i,k)+delp(i,k+1)) qmax = max(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1)) - a4(1,i,k) qmin = a4(1,i,k) - min(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1)) dc(i,k) = sign(min(abs(tmp),qmax,qmin), tmp) df2(i,k) = tmp enddo enddo !----------------------------------------------------------- ! 4th order interpolation of the provisional cell edge value !----------------------------------------------------------- do k=max(3,kmap), km1 do i=i1,i2 c1 = delq(i,k-1)*delp(i,k-1) / d4(i,k) a1 = d4(i,k-1) / (d4(i,k) + delp(i,k-1)) a2 = d4(i,k+1) / (d4(i,k) + delp(i,k)) a4(2,i,k) = a4(1,i,k-1) + c1 + 2./(d4(i,k-1)+d4(i,k+1)) * & ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) - & delp(i,k-1)*a1*dc(i,k ) ) enddo enddo if(km>8 .and. kord>3) call steepz_sak(i1, i2, km, kmap, a4, df2, dc, delq, delp, d4) ! Area preserving cubic with 2nd deriv. = 0 at the boundaries ! Top if ( kmap <= 2 ) then do i=i1,i2 d1 = delp(i,1) d2 = delp(i,2) qm = (d2*a4(1,i,1)+d1*a4(1,i,2)) / (d1+d2) dq = 2.*(a4(1,i,2)-a4(1,i,1)) / (d1+d2) c1 = 4.*(a4(2,i,3)-qm-d2*dq) / ( d2*(2.*d2*d2+d1*(d2+3.*d1)) ) c3 = dq - 0.5*c1*(d2*(5.*d1+d2)-3.*d1**2) a4(2,i,2) = qm - 0.25*c1*d1*d2*(d2+3.*d1) a4(2,i,1) = d1*(2.*c1*d1**2-c3) + a4(2,i,2) dc(i,1) = a4(1,i,1) - a4(2,i,1) ! No over- and undershoot condition cmax = max(a4(1,i,1), a4(1,i,2)) cmin = min(a4(1,i,1), a4(1,i,2)) a4(2,i,2) = max(cmin,a4(2,i,2)) a4(2,i,2) = min(cmax,a4(2,i,2)) enddo endif ! Bottom ! Area preserving cubic with 2nd deriv. = 0 at the surface do i=i1,i2 d1 = delp(i,km) d2 = delp(i,km1) qm = (d2*a4(1,i,km)+d1*a4(1,i,km1)) / (d1+d2) dq = 2.*(a4(1,i,km1)-a4(1,i,km)) / (d1+d2) c1 = (a4(2,i,km1)-qm-d2*dq) / (d2*(2.*d2*d2+d1*(d2+3.*d1))) c3 = dq - 2.0*c1*(d2*(5.*d1+d2)-3.*d1**2) a4(2,i,km) = qm - c1*d1*d2*(d2+3.*d1) a4(3,i,km) = d1*(8.*c1*d1**2-c3) + a4(2,i,km) dc(i,km) = a4(3,i,km) - a4(1,i,km) ! No over- and under-shoot condition cmax = max(a4(1,i,km), a4(1,i,km1)) cmin = min(a4(1,i,km), a4(1,i,km1)) a4(2,i,km) = max(cmin,a4(2,i,km)) a4(2,i,km) = min(cmax,a4(2,i,km)) enddo do k=max(1,kmap),km1 do i=i1,i2 a4(3,i,k) = a4(2,i,k+1) enddo enddo ! Enforce monotonicity of the "slope" within the top layer if ( kmap <= 2 ) then do i=i1,i2 if ( a4(2,i,1) * a4(1,i,1) <= 0. ) then a4(2,i,1) = 0. dc(i,1) = a4(1,i,1) endif if ( dc(i,1) * (a4(2,i,2) - a4(1,i,1)) <= 0. ) then ! Setting DC==0 will force piecewise constant distribution after ! calling kmppm_sak dc(i,1) = 0. endif enddo endif ! Enforce constraint on the "slope" at the surface do i=i1,i2 if( a4(3,i,km) * a4(1,i,km) <= 0. ) then ! a4(3,i,km) = 0. ! dc(i,km) = -a4(1,i,km) dc(i,km) = 0. endif if( dc(i,km) * (a4(1,i,km) - a4(2,i,km)) <= 0. ) then dc(i,km) = 0. endif enddo !----------------------------------------------------------- ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 ) !----------------------------------------------------------- ! Top 2 and bottom 2 layers always use monotonic mapping if ( kmap <= 2 ) then do k=1,2 do i=i1,i2 a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) enddo call kmppm_sak(dc(i1,k), a4(1,i1,k), it, 0) enddo endif if(kord >= 7) then !----------------------- ! Huynh's 2nd constraint !----------------------- do k=max(2,kmap-1), km1 do i=i1,i2 ! Method#1 ! h2(i,k) = delq(i,k) - delq(i,k-1) ! Method#2 ! h2(i,k) = 2.*(dc(i,k+1)/delp(i,k+1) - dc(i,k-1)/delp(i,k-1)) ! & / ( delp(i,k)+0.5*(delp(i,k-1)+delp(i,k+1)) ) ! & * delp(i,k)**2 ! Method#3 h2(i,k) = dc(i,k+1) - dc(i,k-1) enddo enddo if( kord == 7 ) then fac = 1.5 ! original quasi-monotone else fac = 0.125 ! full monotone endif do k=max(3,kmap), km-2 do i=i1,i2 ! Right edges ! qmp = a4(1,i,k) + 2.0*delq(i,k-1) ! lac = a4(1,i,k) + fac*h2(i,k-1) + 0.5*delq(i,k-1) ! pmp = 2.*dc(i,k) qmp = a4(1,i,k) + pmp lac = a4(1,i,k) + fac*h2(i,k-1) + dc(i,k) qmin = min(a4(1,i,k), qmp, lac) qmax = max(a4(1,i,k), qmp, lac) a4(3,i,k) = min(max(a4(3,i,k), qmin), qmax) ! Left edges ! qmp = a4(1,i,k) - 2.0*delq(i,k) ! lac = a4(1,i,k) + fac*h2(i,k+1) - 0.5*delq(i,k) ! qmp = a4(1,i,k) - pmp lac = a4(1,i,k) + fac*h2(i,k+1) - dc(i,k) qmin = min(a4(1,i,k), qmp, lac) qmax = max(a4(1,i,k), qmp, lac) a4(2,i,k) = min(max(a4(2,i,k), qmin), qmax) ! Recompute A6 a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) enddo ! Additional constraint to ensure positivity when kord=7 if (iv == 0 .and. kord == 7) then call kmppm_sak(dc(i1,k), a4(1,i1,k), it, 2) endif enddo else lmt = kord - 3 lmt = max(0, lmt) if (iv == 0) lmt = min(2, lmt) do k=max(3,kmap), km-2 if( kord /= 4) then do i=i1,i2 a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) enddo endif call kmppm_sak(dc(i1,k), a4(1,i1,k), it, lmt) enddo endif do k=km1,km do i=i1,i2 a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) enddo call kmppm_sak(dc(i1,k), a4(1,i1,k), it, 0) enddo !EOC end subroutine ppm2m_sak !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! !ROUTINE: kmppm_sak --- Perform piecewise parabolic method in vertical ! ! !INTERFACE: subroutine kmppm_sak(dm, a4, itot, lmt) implicit none ! !INPUT PARAMETERS: real, intent(in):: dm(*) ! the linear slope integer, intent(in) :: itot ! Total Longitudes integer, intent(in) :: lmt ! 0: Standard PPM constraint ! 1: Improved full monotonicity constraint (Lin) ! 2: Positive definite constraint ! 3: do nothing (return immediately) ! !INPUT/OUTPUT PARAMETERS: real, intent(inout) :: a4(4,*) ! PPM array ! AA <-- a4(1,i) ! AL <-- a4(2,i) ! AR <-- a4(3,i) ! A6 <-- a4(4,i) ! !DESCRIPTION: ! ! !REVISION HISTORY: ! 00.04.24 Lin Last modification ! 01.03.26 Sawyer Added ProTeX documentation ! 02.04.04 Sawyer Incorporated newest FVGCM version ! !EOP !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: real, parameter:: r12 = 1./12. real qmp real da1, da2, a6da real fmin integer i ! Developer: S.-J. Lin, NASA-GSFC ! Last modified: Apr 24, 2000 if ( lmt == 3 ) return if(lmt == 0) then ! Standard PPM constraint do i=1,itot if(dm(i) == 0.) then a4(2,i) = a4(1,i) a4(3,i) = a4(1,i) a4(4,i) = 0. else da1 = a4(3,i) - a4(2,i) da2 = da1**2 a6da = a4(4,i)*da1 if(a6da < -da2) then a4(4,i) = 3.*(a4(2,i)-a4(1,i)) a4(3,i) = a4(2,i) - a4(4,i) elseif(a6da > da2) then a4(4,i) = 3.*(a4(3,i)-a4(1,i)) a4(2,i) = a4(3,i) - a4(4,i) endif endif enddo elseif (lmt == 1) then ! Improved full monotonicity constraint (Lin 2003) ! Note: no need to provide first guess of A6 <-- a4(4,i) do i=1, itot qmp = 2.*dm(i) a4(2,i) = a4(1,i)-sign(min(abs(qmp),abs(a4(2,i)-a4(1,i))), qmp) a4(3,i) = a4(1,i)+sign(min(abs(qmp),abs(a4(3,i)-a4(1,i))), qmp) a4(4,i) = 3.*( 2.*a4(1,i) - (a4(2,i)+a4(3,i)) ) enddo elseif (lmt == 2) then ! Positive definite constraint do i=1,itot if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then fmin = a4(1,i)+0.25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12 if( fmin < 0. ) then if(a4(1,i) a4(2,i)) then a4(4,i) = 3.*(a4(2,i)-a4(1,i)) a4(3,i) = a4(2,i) - a4(4,i) else a4(4,i) = 3.*(a4(3,i)-a4(1,i)) a4(2,i) = a4(3,i) - a4(4,i) endif endif endif enddo endif !EOC end subroutine kmppm_sak !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! !ROUTINE: steepz_sak --- Calculate attributes for PPM ! ! !INTERFACE: subroutine steepz_sak(i1, i2, km, kmap, a4, df2, dm, dq, dp, d4) implicit none ! !INPUT PARAMETERS: integer, intent(in) :: km ! Total levels integer, intent(in) :: kmap ! integer, intent(in) :: i1 ! Starting longitude integer, intent(in) :: i2 ! Finishing longitude real, intent(in) :: dp(i1:i2,km) ! grid size real, intent(in) :: dq(i1:i2,km) ! backward diff of q real, intent(in) :: d4(i1:i2,km) ! backward sum: dp(k)+ dp(k-1) real, intent(in) :: df2(i1:i2,km) ! first guess mismatch real, intent(in) :: dm(i1:i2,km) ! monotonic mismatch ! !INPUT/OUTPUT PARAMETERS: real, intent(inout) :: a4(4,i1:i2,km) ! first guess/steepened ! ! !DESCRIPTION: ! This is complicated stuff related to the Piecewise Parabolic Method ! and I need to read the Collela/Woodward paper before documenting ! thoroughly. ! ! !REVISION HISTORY: ! ??.??.?? Lin? Creation ! 01.03.26 Sawyer Added ProTeX documentation ! !EOP !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: integer i, k real alfa(i1:i2,km) real f(i1:i2,km) real rat(i1:i2,km) real dg2 ! Compute ratio of dq/dp do k=max(2,kmap-1),km do i=i1,i2 rat(i,k) = dq(i,k-1) / d4(i,k) enddo enddo ! Compute F do k=max(2,kmap-1),km-1 do i=i1,i2 f(i,k) = (rat(i,k+1) - rat(i,k)) & / ( dp(i,k-1)+dp(i,k)+dp(i,k+1) ) enddo enddo do k=max(3,kmap),km-2 do i=i1,i2 if(f(i,k+1)*f(i,k-1)<0. .and. df2(i,k)/=0.) then dg2 = (f(i,k+1)-f(i,k-1))*((dp(i,k+1)-dp(i,k-1))**2 & + d4(i,k)*d4(i,k+1) ) alfa(i,k) = max(0., min(0.5, -0.1875*dg2/df2(i,k))) else alfa(i,k) = 0. endif enddo enddo do k=max(4,kmap+1),km-2 do i=i1,i2 a4(2,i,k) = (1.-alfa(i,k-1)-alfa(i,k)) * a4(2,i,k) + & alfa(i,k-1)*(a4(1,i,k)-dm(i,k)) + & alfa(i,k)*(a4(1,i,k-1)+dm(i,k-1)) enddo enddo !EOC end subroutine steepz_sak !----------------------------------------------------------------------- subroutine cloud_clear_xfer ( cloud_generator_on,tmp3, qa_mean, qa_mean_lst, a_clr, a_cld, clr, cld) logical, intent(in) :: cloud_generator_on real, dimension(:,:), intent(in) :: tmp3, qa_mean, qa_mean_lst real, dimension(:,:), intent(inout) :: a_clr, a_cld, clr, cld real :: cld2clr, clr2cld, prec_cld2clr, prec_clr2cld, tmp1, tmp2 integer :: k,i,kdim,idim idim = size(tmp3,1) kdim = size(tmp3,2) do k=1,kdim do i=1,idim !------------------------------- !compute cloud to clear transfer if (overlap .eq. 1) & cld2clr= min(a_cld(i,k),max(0.,a_cld(i,k) - qa_mean(i,k)) ) if (overlap .eq. 2) & cld2clr= min(a_cld(i,k),max(0.,a_cld(i,k)*(1.-qa_mean(i,k)))) if (cloud_generator_on) then tmp1 = min(a_cld(i,k),max(0.,a_cld(i,k) - qa_mean(i,k)) ) tmp2 = min(a_cld(i,k),max(0.,a_cld(i,k)*(1.-qa_mean(i,k)))) cld2clr=min(a_cld(i,k),max(0.,tmp3(i,k)*tmp1+(1.-tmp3(i,k))*tmp2)) end if !------------------------------- !compute clear to cloud transfer if (overlap .eq. 1) & clr2cld = min(max(qa_mean(i,k)-qa_mean_lst(i,k),0.),a_clr(i,k)) if (overlap .eq. 2) & clr2cld = min(max( a_clr(i,k)*qa_mean(i,k),0.),a_clr(i,k)) if (cloud_generator_on) then tmp1 = min(max(qa_mean(i,k)-qa_mean_lst(i,k),0.),a_clr(i,k)) tmp2 = min(max( a_clr(i,k)*qa_mean(i,k),0.),a_clr(i,k)) clr2cld=min(a_clr(i,k),max(0.,tmp3(i,k)*tmp1+(1.-tmp3(i,k))*tmp2)) end if !--------------------------------- !calculate precipitation transfers prec_cld2clr = cld(i,k)*(cld2clr/max(a_cld(i,k),qmin)) prec_clr2cld = clr(i,k)*(clr2cld/max(a_clr(i,k),qmin)) !---------------- !add in transfers a_clr(i,k) = a_clr(i,k) + cld2clr - clr2cld a_cld(i,k) = a_cld(i,k) - cld2clr + clr2cld clr(i,k) = clr(i,k) + prec_cld2clr - prec_clr2cld cld(i,k) = cld(i,k) - prec_cld2clr + prec_clr2cld enddo enddo end subroutine cloud_clear_xfer ! !####################################################################### !####################################################################### end module strat_cloud_legacy_mod