Index: /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_globals.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_globals.F90 (revision 3951) +++ /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_globals.F90 (revision 3952) @@ -30,13 +30,17 @@ ! - Vertical structure part ! - ! The module also contains thirteen methods: + ! The module also contains twenty methods: ! - mm_global_init_0 ! - mm_global_init_1 ! - mm_column_init ! - mm_aerosols_init + ! - mm_clouds_init ! - mm_alpha_s, mm_alpha_f - ! - mm_effg ! - mm_set_moments_thresholds ! - mm_get_rcs, mm_get_rcf + ! - mm_set_moments_cld_thresholds + ! - cldprop_sc, cldprop_ve + ! - read_esp + ! - mm_effg ! - mm_dump_parameters ! - check_r1, check_l1, check_i1, check_s1 @@ -58,5 +62,5 @@ PUBLIC - PRIVATE :: check_r1,check_i1,check_l1,check_s1,read_esp + PRIVATE :: cldprop_sc,cldprop_ve,read_esp,check_r1,check_i1,check_l1,check_s1 ! ~~~~~~~~~~~~~~~~~~~ @@ -71,5 +75,5 @@ PROTECTED :: mm_ini,mm_ini_col,mm_ini_aer ! Model parameters (mm_global_init) - PROTECTED :: mm_dt,mm_rhoaer,mm_df,mm_rm,mm_haze_prod_pCH4,mm_p_prod,mm_rc_prod,mm_tx_prod,mm_rpla,mm_g0,mm_rb2ra + PROTECTED :: mm_dt,mm_rhoaer,mm_df,mm_rm,mm_call_CH4hazeprod,mm_p_prod,mm_rc_prod,mm_tx_prod,mm_rpla,mm_g0,mm_rb2ra ! Atmospheric vertical structure (mm_column_init) PROTECTED :: mm_nla,mm_nle,mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay @@ -84,13 +88,15 @@ !~~~~~~~~~~~~~ ! Initialization control flags (cannot be updated) - ! PROTECTED :: mm_ini_cld + PROTECTED :: mm_ini_cld ! Condensible species parameters (mm_global_init) PROTECTED :: mm_nesp, mm_spcname, mm_xESPS + ! Condensible species parameters (mm_clouds_init) + PROTECTED :: mm_gas ! Moments parameters (mm_clouds_init) - ! PROTECTED :: mm_m0ccn, mm_m3ccn, mm_m3ice + PROTECTED :: mm_m0ccn, mm_m3ccn, mm_m3ice ! Moments parameters (derived, are updated with moments parameters) - ! PROTECTED :: mm_drad, mm_drho + PROTECTED :: mm_drad, mm_drho ! Thresholds parameters - ! PROTECTED :: mm_m0ccn_min, mm_m3cld_min + PROTECTED :: mm_m0ccn_min, mm_m3ccn_min, mm_m3cld_min, mm_drad_min, mm_drad_max ! ~~~~~~~~~~~~~ @@ -100,19 +106,19 @@ !~~~~~~~~~~~~ ! Enable/Disable haze production. - LOGICAL, SAVE :: mm_w_haze_prod = .true. + LOGICAL, SAVE :: mm_call_hazeprod = .true. ! Enable/Disable haze production from CH4 photolysis. - LOGICAL, SAVE :: mm_haze_prod_pCH4 = .true. + LOGICAL, SAVE :: mm_call_CH4hazeprod = .true. ! Enable/Disable haze sedimentation. - LOGICAL, SAVE :: mm_w_haze_sed = .true. + LOGICAL, SAVE :: mm_call_hazesed = .true. ! Enable/Disable haze coagulation. - LOGICAL, SAVE :: mm_w_haze_coag = .true. + LOGICAL, SAVE :: mm_call_hazecoag = .true. ! Force all aerosols moments to fall at M0 settling velocity. - LOGICAL, SAVE :: mm_wsed_m0 = .false. + LOGICAL, SAVE :: mm_wsed_m0 = .false. ! Force all aerosols moments to fall at M3 settling velocity. - LOGICAL, SAVE :: mm_wsed_m3 = .false. + LOGICAL, SAVE :: mm_wsed_m3 = .false. ! Enable/Disable spherical probability transfert. - LOGICAL, SAVE :: mm_w_ps2s = .true. + LOGICAL, SAVE :: mm_call_ps2s = .true. ! Enable/Disable aerosol electric charge correction. - LOGICAL, SAVE :: mm_w_qe = .true. + LOGICAL, SAVE :: mm_call_qe = .true. ! Cloud model: @@ -134,4 +140,6 @@ ! Initialization control flag [[mm_globals(module):mm_aerosols_init(function)]]. LOGICAL, PUBLIC, SAVE :: mm_ini_aer = .false. + ! Initialization control flag [[mm_globals(module):mm_clouds_init(function)]]. + LOGICAL, PUBLIC, SAVE :: mm_ini_cld = .false. ! ~~~~~~~~~~~~~~~~~~~~~~~ @@ -197,11 +205,13 @@ !~~~~~~~~~~~~~ ! Total number of cloud condensation nuclei minimum threshold. - ! REAL(kind=mm_wp), SAVE :: mm_m0ccn_min = 1.e-8_mm_wp + REAL(kind=mm_wp), SAVE :: mm_m0ccn_min = 1.e-8_mm_wp + ! Total volume of cloud condensation nuclei minimum threshold. + REAL(kind=mm_wp), SAVE :: mm_m3ccn_min = 1.e-35_mm_wp ! Total volume of cloud drop minimum threshold. - ! REAL(kind=mm_wp), SAVE :: mm_m3cld_min = 1.e-35_mm_wp + REAL(kind=mm_wp), SAVE :: mm_m3cld_min = 1.e-35_mm_wp ! Characteristic cloud drop radius minimum threshold. - ! REAL(kind=mm_wp), SAVE :: mm_drad_min = 1.e-9_mm_wp + REAL(kind=mm_wp), SAVE :: mm_drad_min = 1.e-9_mm_wp ! Characteristic cloud drop radius Maximum threshold. - ! REAL(kind=mm_wp), SAVE :: mm_drad_max = 1.e-2_mm_wp + REAL(kind=mm_wp), SAVE :: mm_drad_max = 1.e-3_mm_wp ! Planet radius (m) and gravity acceleration (m.s-2). @@ -354,4 +364,6 @@ REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlay + ! Haze model: + ! ~~~~~~~~~~~ ! Spherical mode - 0th order moment (m-3). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_s @@ -392,13 +404,57 @@ REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_f_flux + ! Cloud model: + ! ~~~~~~~~~~~~ + ! Cloud condensation nuclei - 0th order moment (m-3). + REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0ccn + ! Cloud condensation nuclei - 3rd order moment (m3.m-3). + REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3ccn + ! Mean drop radius (m). + REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drad + ! Mean Drop density (kg.m-3). + REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drho + + ! Ice components - 3rd order moments (m3.m-3). + ! It is a 2D array with the vertical layers in first dimension, and + ! the number of ice components in the second. + REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_m3ice + ! Condensible species molar fraction (mol.mol-1). + ! It is a 2D array with the vertical layers in first dimension, and + ! the number of condensible species in the second. + REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gas + + ! CCN precipitation (kg.m-2). + REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp + ! Ice components precipitation (kg.m-2). + REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ice_prec + + ! CCN (and ices) settling velocity (m.s-1). + REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_cld_vsed + + ! CCN mass fluxes (kg.m-2.s-1). + REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_flux + ! Ice components sedimentation fluxes (kg.m-2.s-1). + REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_ice_fluxes + + ! Condensible species saturation ratio. + REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gas_sat + ! Condensible components nucleation rates (m-2.s-1). + REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_nrate + ! Condensible components growth rates (m2.s-1). + REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_grate + + ! --- OPENMP --------------- ! All variables related to column computations should be private to each thread - !$OMP THREADPRIVATE(mm_ini_col,mm_ini_aer) + !$OMP THREADPRIVATE(mm_ini_col,mm_ini_aer,mm_ini_cld) !$OMP THREADPRIVATE(mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay) !$OMP THREADPRIVATE(mm_m0aer_s,mm_m3aer_s,mm_m0aer_f,mm_m3aer_f) - !$OMP THREADPRIVATE(mm_rcs,mm_rcf) - !$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed) - !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux) + !$OMP THREADPRIVATE(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_gas) + !$OMP THREADPRIVATE(mm_rcs,mm_rcf,mm_drad,mm_drho) + !$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed,mm_cld_vsed) + !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_prec,mm_ice_prec,mm_ccn_flux,mm_ice_fluxes,mm_gas_sat,mm_nrate,mm_grate) !$OMP THREADPRIVATE(mm_m0as_min,mm_m3as_min,mm_rcs_min,mm_m0af_min,mm_m3af_min,mm_rcf_min) + !$OMP THREADPRIVATE(mm_m0ccn_min,mm_m3ccn_min,mm_m3cld_min,mm_drad_min,mm_drad_max) !$OMP THREADPRIVATE(mm_nla,mm_nle) + ! -------------------------- ! Interface to global initialization. @@ -417,4 +473,12 @@ END INTERFACE mm_check_opt + ! Interface to cloud properties methods. + ! The method computes clouds properties (mean drop radius and denstity) from their afferent + ! moments. It is overloaded to compute properties at a single level or over all the vertical + ! atmospheric structure. + INTERFACE mm_cloud_properties + MODULE PROCEDURE cldprop_sc,cldprop_ve + END INTERFACE mm_cloud_properties + CONTAINS @@ -425,10 +489,10 @@ !============================================================================ - FUNCTION mm_global_init_0(dt,df,rm,rho_aer,haze_prod_pCH4,p_prod,tx_prod,rc_prod,& - rplanet,g0,air_rad,air_mmol, & - coag_interactions,w_haze_prod,w_haze_sed,w_haze_coag, & - force_wsed_to_m0,force_wsed_to_m3, & - m0as_min,rcs_min,m0af_min,rcf_min, & - clouds,spcfile,debug) RESULT(err) + FUNCTION mm_global_init_0(dt,df,rm,rho_aer,call_CH4hazeprod,p_prod,tx_prod,rc_prod, & + rplanet,g0,air_rad,air_mmol, & + coag_interactions,call_hazeprod,call_hazesed,call_hazecoag,& + force_wsed_to_m0,force_wsed_to_m3, & + m0as_min,rcs_min,m0af_min,rcf_min, & + clouds,spcfile,m0ccn_min,drad_min,debug) RESULT(err) !! Initialize global parameters of the model. !! @@ -449,5 +513,5 @@ REAL(kind=mm_wp), INTENT(in) :: rho_aer ! Enable/Disable production from CH4 photolysis. - LOGICAL, INTENT(in) :: haze_prod_pCH4 + LOGICAL, INTENT(in) :: call_CH4hazeprod ! Aerosol production pressure level (Pa). REAL(kind=mm_wp), INTENT(in) :: p_prod @@ -468,9 +532,9 @@ INTEGER, INTENT(in) :: coag_interactions ! Haze production process control flag. - LOGICAL, INTENT(in) :: w_haze_prod + LOGICAL, INTENT(in) :: call_hazeprod ! Haze sedimentation process control flag. - LOGICAL, INTENT(in) :: w_haze_sed + LOGICAL, INTENT(in) :: call_hazesed ! Haze coagulation process control flag. - LOGICAL, INTENT(in) :: w_haze_coag + LOGICAL, INTENT(in) :: call_hazecoag ! Force __all__ aerosols moments to fall at M0 settling velocity. LOGICAL, INTENT(in) :: force_wsed_to_m0 @@ -491,4 +555,8 @@ ! Clouds microphysics condensible species properties file. CHARACTER(len=*), INTENT(in) :: spcfile + ! Minimum threshold for M0 of cloud condensation nuceli (m-3). + REAL(kind=mm_wp), INTENT(in) :: m0ccn_min + ! Minimum threshold for the cloud drop radius (m). + REAL(kind=mm_wp), INTENT(in) :: drad_min ! Debug mode control flag. @@ -512,17 +580,17 @@ ! Free parameters: - mm_dt = dt - mm_df = df - mm_rm = rm - mm_rhoaer = rho_aer - mm_haze_prod_pCH4 = haze_prod_pCH4 - mm_p_prod = p_prod - mm_tx_prod = tx_prod - mm_rc_prod = rc_prod - mm_rpla = rplanet - mm_g0 = g0 - mm_air_rad = air_rad - mm_air_mmol = air_mmol - mm_air_mmas = air_mmol / mm_navo + mm_dt = dt + mm_df = df + mm_rm = rm + mm_rhoaer = rho_aer + mm_call_CH4hazeprod = call_CH4hazeprod + mm_p_prod = p_prod + mm_tx_prod = tx_prod + mm_rc_prod = rc_prod + mm_rpla = rplanet + mm_g0 = g0 + mm_air_rad = air_rad + mm_air_mmol = air_mmol + mm_air_mmas = air_mmol / mm_navo ! Microphysical processes - Haze: @@ -532,13 +600,13 @@ RETURN ENDIF - mm_w_haze_prod = w_haze_prod - mm_w_haze_sed = w_haze_sed - mm_w_haze_coag = w_haze_coag - mm_wsed_m0 = force_wsed_to_m0 - mm_wsed_m3 = force_wsed_to_m3 + mm_call_hazeprod = call_hazeprod + mm_call_hazesed = call_hazesed + mm_call_hazecoag = call_hazecoag + mm_wsed_m0 = force_wsed_to_m0 + mm_wsed_m3 = force_wsed_to_m3 ! Moment threshold flags: mm_m0as_min = MAX(0._mm_wp,m0as_min) - mm_rcs_min = MAX(1.e-10_mm_wp,rcs_min) + mm_rcs_min = MAX(1.e-9_mm_wp,rcs_min) mm_m0af_min = MAX(0._mm_wp,m0af_min) mm_rcf_min = MAX(mm_rm,rcf_min) @@ -591,4 +659,12 @@ ENDIF ENDDO + + ! Moment threshold flags: + mm_m0ccn_min = MAX(0._mm_wp,m0ccn_min) + mm_drad_min = MAX(1.e-9_mm_wp,drad_min) + + ! Computes M3 thresholds: + mm_m3ccn_min = mm_m0ccn_min*mm_alpha_s(3._mm_wp) * mm_rcs_min**3._mm_wp + mm_m3cld_min = mm_m0ccn_min * mm_drad_min**3._mm_wp ENDIF ! end of mm_call_clouds @@ -645,5 +721,5 @@ err = mm_check_opt(cfg_get_value(cfg,"rho_aer",mm_rhoaer),mm_rhoaer,wlog=mm_log) IF (err/=0) RETURN - err = mm_check_opt(cfg_get_value(cfg,"call_haze_prod",mm_haze_prod_pCH4),mm_haze_prod_pCH4,wlog=mm_log) + err = mm_check_opt(cfg_get_value(cfg,"call_haze_prod",mm_call_CH4hazeprod),mm_call_CH4hazeprod,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"p_prod",mm_p_prod),mm_p_prod,wlog=mm_log) @@ -671,9 +747,9 @@ RETURN ENDIF - err = mm_check_opt(cfg_get_value(cfg,"haze_production",mm_w_haze_prod),mm_w_haze_prod,wlog=mm_log) - IF (err/=0) RETURN - err = mm_check_opt(cfg_get_value(cfg,"haze_sedimentation",mm_w_haze_sed),mm_w_haze_sed,wlog=mm_log) - IF (err/=0) RETURN - err = mm_check_opt(cfg_get_value(cfg,"haze_coagulation",mm_w_haze_coag),mm_w_haze_coag,wlog=mm_log) + err = mm_check_opt(cfg_get_value(cfg,"haze_production",mm_call_hazeprod),mm_call_hazeprod,wlog=mm_log) + IF (err/=0) RETURN + err = mm_check_opt(cfg_get_value(cfg,"haze_sedimentation",mm_call_hazesed),mm_call_hazesed,wlog=mm_log) + IF (err/=0) RETURN + err = mm_check_opt(cfg_get_value(cfg,"haze_coagulation",mm_call_hazecoag),mm_call_hazecoag,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"wsed_m0",mm_wsed_m0),mm_wsed_m0,wlog=mm_log) @@ -936,8 +1012,8 @@ ! Initialization of aerosol tracers: ! @note: mm_dzlev is already from top to ground - mm_m0aer_s = m0aer_s(mm_nla:1:-1)/mm_dzlev(:) - mm_m3aer_s = m3aer_s(mm_nla:1:-1)/mm_dzlev(:) - mm_m0aer_f = m0aer_f(mm_nla:1:-1)/mm_dzlev(:) - mm_m3aer_f = m3aer_f(mm_nla:1:-1)/mm_dzlev(:) + mm_m0aer_s = m0aer_s(mm_nla:1:-1) / mm_dzlev(:) + mm_m3aer_s = m3aer_s(mm_nla:1:-1) / mm_dzlev(:) + mm_m0aer_f = m0aer_f(mm_nla:1:-1) / mm_dzlev(:) + mm_m3aer_f = m3aer_f(mm_nla:1:-1) / mm_dzlev(:) ! Setup threshold (useful for preventing bugs): @@ -962,4 +1038,89 @@ END FUNCTION mm_aerosols_init + + FUNCTION mm_clouds_init(m0ccn,m3ccn,m3ice,gas) RESULT(err) + !! Initialize clouds tracers vertical grid. + !! + !! The subroutine initializes cloud microphysics tracers columns. It allocates variables if + !! required and stores input vectors in reversed order. It also computes the mean drop radius + !! and density and allocates diagnostic vectors. + !! + !! @note + !! All the input arguments should be defined from __GROUND__ to __TOP__. + !! + !! @warning + !! [[mm_global_init(interface)]] and [[mm_column_init(function)]] must have been called at least once before this method is called. + !! Moreover, this method should be after each call of [[mm_column_init(function)]] to reflect changes in the + !! vertical atmospheric structure. + !! + REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0ccn ! 0th order moment of the CCN distribution (m-2). + REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3ccn ! 3rd order moment of the CCN distribution (m3.m-2). + REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: m3ice ! 3rd order moments of the ice components (m3.m-2). + REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: gas ! Condensible species gas molar fraction (mol.mol-1). + + INTEGER :: i + TYPE(error) :: err ! Error status of the function. + + err = noerror + IF (.NOT.mm_ini) THEN + err = error("Global initialization not done yet",-8) + RETURN + ENDIF + + IF (.NOT.mm_call_clouds) THEN + IF (mm_debug) WRITE(*,'(a)') "[MM_DEBUG - mm_clouds_init]: Cloud microphysic is not enabled..." + RETURN + ENDIF + + ! Allocate variable if required: + IF (.NOT.ALLOCATED(mm_m0ccn)) ALLOCATE(mm_m0ccn(mm_nla)) + IF (.NOT.ALLOCATED(mm_m3ccn)) ALLOCATE(mm_m3ccn(mm_nla)) + IF (.NOT.ALLOCATED(mm_m3ice)) ALLOCATE(mm_m3ice(mm_nla,mm_nesp)) + IF (.NOT.ALLOCATED(mm_gas)) ALLOCATE(mm_gas(mm_nla,mm_nesp)) + IF (.NOT.ALLOCATED(mm_drad)) ALLOCATE(mm_drad(mm_nla)) + IF (.NOT.ALLOCATED(mm_drho)) ALLOCATE(mm_drho(mm_nla)) + + ! Allocate memory for diagnostics: + mm_ccn_prec = 0._mm_wp + IF (.NOT.ALLOCATED(mm_ice_prec)) THEN + ALLOCATE(mm_ice_prec(mm_nesp)) ; mm_ice_prec(:) = 0._mm_wp + ENDIF + IF (.NOT.ALLOCATED(mm_cld_vsed)) THEN + ALLOCATE(mm_cld_vsed(mm_nla)) ; mm_cld_vsed(:) = 0._mm_wp + ENDIF + IF (.NOT.ALLOCATED(mm_ccn_flux)) THEN + ALLOCATE(mm_ccn_flux(mm_nla)) ; mm_ccn_flux(:) = 0._mm_wp + ENDIF + IF (.NOT.ALLOCATED(mm_ice_fluxes)) THEN + ALLOCATE(mm_ice_fluxes(mm_nla,mm_nesp)) ; mm_ice_fluxes(:,:) = 0._mm_wp + ENDIF + IF (.NOT.ALLOCATED(mm_gas_sat)) THEN + ALLOCATE(mm_gas_sat(mm_nla,mm_nesp)) ; mm_gas_sat(:,:) = 0._mm_wp + ENDIF + IF (.NOT.ALLOCATED(mm_nrate)) THEN + ALLOCATE(mm_nrate(mm_nla,mm_nesp)) ; mm_nrate(:,:) = 0._mm_wp + ENDIF + IF (.NOT.ALLOCATED(mm_grate)) THEN + ALLOCATE(mm_grate(mm_nla,mm_nesp)) ; mm_grate(:,:) = 0._mm_wp + ENDIF + + ! /!\ mm_dzlev already from top to ground + mm_m0ccn = m0ccn(mm_nla:1:-1) / mm_dzlev(:) + mm_m3ccn = m3ccn(mm_nla:1:-1) / mm_dzlev(:) + DO i = 1, mm_nesp + mm_m3ice(:,i) = m3ice(mm_nla:1:-1,i) / mm_dzlev(:) + mm_gas(:,i) = gas(mm_nla:1:-1,i) + ENDDO + + ! Setup threshold: + call mm_set_moments_cld_thresholds() + + ! Drop mean radius and density: + call mm_cloud_properties(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_drad,mm_drho) + + ! End of initialization: + mm_ini_cld = .true. + + END FUNCTION mm_clouds_init @@ -1063,4 +1224,95 @@ ! CLOUD RELATED METHODS !============================================================================ + + SUBROUTINE mm_set_moments_cld_thresholds() + !! Apply minimum threshold for the cloud drop moments. + !! The method resets moments values to zero if their current value is below the minimum threholds. + !! + INTEGER :: i, j + REAL(kind=mm_wp) :: m3cld = 0._mm_wp + + DO i = 1, mm_nla + m3cld = mm_m3ccn(i) + DO j = 1, mm_nesp + m3cld = m3cld + mm_m3ice(i,j) + ENDDO + + IF ((mm_m0ccn(i) < mm_m0ccn_min) .OR. (mm_m3ccn(i) < mm_m3ccn_min) .OR. (m3cld < mm_m3cld_min)) THEN + mm_m0ccn(i) = 0._mm_wp + mm_m3ccn(i) = 0._mm_wp + DO j = 1, mm_nesp + mm_m3ice(i,j) = 0._mm_wp + ENDDO + ENDIF + ENDDO + END SUBROUTINE mm_set_moments_cld_thresholds + + + SUBROUTINE cldprop_sc(m0ccn,m3ccn,m3ice,drad,drho) + !! Get cloud drop properties (scalar). + !! The method computes the mean radius and mean density of cloud drops. + !! + !! @warning + !! If __drad__ is greater than __drmax__ it is automatically set to __drmax__, but computation of + !! __drho__ remains unmodified. So __drho__ is not correct in that case! + !! + REAL(kind=mm_wp), INTENT(in) :: m0ccn ! 0th order moment of the ccn + REAL(kind=mm_wp), INTENT(in) :: m3ccn ! 3rd order moment of the ccn + REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ice ! 3rd order moments of each ice component + REAL(kind=mm_wp), INTENT(out) :: drad ! Output mean drop radius + REAL(kind=mm_wp), INTENT(out), OPTIONAL :: drho ! Optional output mean drop density + + REAL(kind=mm_wp) :: Ntot, Vtot, Wtot + REAL(kind=mm_wp), PARAMETER :: athird = 1._mm_wp / 3._mm_wp + REAL(kind=mm_wp), PARAMETER :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp + + ! Set to zero : + drad = 0._mm_wp + IF (PRESENT(drho)) drho = 0._mm_wp + + ! Initialization : + Ntot = m0ccn + Vtot = m3ccn + SUM(m3ice) + Wtot = m3ccn*mm_rhoaer + SUM(m3ice*mm_xESPS(:)%rho_s) + + IF (Ntot <= mm_m0ccn_min .OR. Vtot <= mm_m3cld_min) THEN + drad = mm_drad_min + IF (PRESENT(drho)) drho = mm_rhoaer + ELSE + drad = (Vtot / Ntot)**athird + drad = MAX(MIN(drad,mm_drad_max),mm_drad_min) + IF (PRESENT(drho)) drho = Wtot / Vtot + ENDIF + + RETURN + END SUBROUTINE cldprop_sc + + SUBROUTINE cldprop_ve(m0ccn,m3ccn,m3ice,drad,drho) + !! Get cloud drop properties (vector). + !! + !! The method performs the same computations than [[cldprop_sc(subroutine)]] + !! but for the entire vertical atmospheric structure. + !! + REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m0ccn ! 0th order moment of the ccn. + REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ccn ! 3rd order moment of the ccn. + REAL(kind=mm_wp), INTENT(in), DIMENSION(:,:) :: m3ice ! 3rd order moments of each ice component. + REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: drad ! Output mean drop radius. + REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: drho ! Optional output mean drop density. + + INTEGER :: i + + IF (PRESENT(drho)) THEN + DO i = 1, SIZE(m0ccn) + call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i),drho(i)) + ENDDO + ELSE + DO i = 1, SIZE(m0ccn) + call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i)) + ENDDO + ENDIF + + RETURN + END SUBROUTINE cldprop_ve + FUNCTION read_esp(parser,sec,pp) RESULT (err) @@ -1127,8 +1379,8 @@ ! Haze production: - WRITE(*,'(a,L2)') "mm_w_haze_prod : ", mm_w_haze_prod + WRITE(*,'(a,L2)') "mm_call_hazeprod : ", mm_call_hazeprod WRITE(*,'(a,ES14.7)') "mm_rc_prod (m) : ", mm_rc_prod - WRITE(*,'(a,L2)') "mm_haze_prod_pCH4 : ", mm_haze_prod_pCH4 - IF (.NOT. mm_haze_prod_pCH4) THEN + WRITE(*,'(a,L2)') "mm_call_CH4hazeprod : ", mm_call_CH4hazeprod + IF (.NOT. mm_call_CH4hazeprod) THEN WRITE(*,'(a,ES14.7)') " --> mm_p_prod (Pa) : ", mm_p_prod WRITE(*,'(a,ES14.7)') " --> mm_tx_prod (kg.m-2.s-1): ", mm_tx_prod @@ -1137,12 +1389,12 @@ ! Haze coagulation: - WRITE(*,'(a,L2)') "mm_w_haze_coag : ", mm_w_haze_coag - IF (mm_w_haze_coag) THEN + WRITE(*,'(a,L2)') "mm_call_hazecoag : ", mm_call_hazecoag + IF (mm_call_hazecoag) THEN WRITE(*,'(a,I2.2)') " --> mm_coag_interactions : ", mm_coag_choice ENDIF ! Haze sedimentation: - WRITE(*,'(a,L2)') "mm_w_haze_sed : ", mm_w_haze_sed - IF (mm_w_haze_sed) THEN + WRITE(*,'(a,L2)') "mm_call_hazesed : ", mm_call_hazesed + IF (mm_call_hazesed) THEN WRITE(*,'(a,L2)') " --> mm_wsed_m0 : ", mm_wsed_m0 WRITE(*,'(a,L2)') " --> mm_wsed_m3 : ", mm_wsed_m3 @@ -1157,4 +1409,14 @@ WRITE(*,'(a,ES14.7)') " mm_m0af_min (m-3) : ", mm_m0af_min WRITE(*,'(a,ES14.7)') " mm_rcf_min (m) : ", mm_rcf_min + WRITE(*,'(a)') "---------------------------------------" + + ! Clouds related: + WRITE(*,'(a,L2)') "mm_call_clouds : ", mm_call_clouds + IF (mm_call_clouds) THEN + WRITE(*,'(a)') " Thresholds clouds drop" + WRITE(*,'(a,ES14.7)') " --> mm_m0ccn_min : ", mm_m0ccn_min + WRITE(*,'(a,ES14.7)') " --> mm_drad_min : ", mm_drad_min + WRITE(*,'(a,ES14.7)') " --> mm_drad_max : ", mm_drad_max + ENDIF WRITE(*,'(a)') "---------------------------------------" Index: /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_haze.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_haze.F90 (revision 3951) +++ /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_haze.F90 (revision 3952) @@ -99,10 +99,10 @@ ! Coagulation - IF (mm_w_haze_coag) THEN + IF (mm_call_hazecoag) THEN call mm_haze_coagulation(dm0a_s,dm3a_s,dm0a_f,dm3a_f) ENDIF ! Sedimentation - IF (mm_w_haze_sed) THEN + IF (mm_call_hazesed) THEN call mm_haze_sedimentation(zdm0as,zdm3as,zdm0af,zdm3af) @@ -119,7 +119,7 @@ ! Production - IF (mm_w_haze_prod) THEN + IF (mm_call_hazeprod) THEN ! Production by photolysis of CH4 - IF (mm_haze_prod_pCH4) THEN + IF (mm_call_CH4hazeprod) THEN dm0a_s = dm0a_s + (m3a_s_prod / (mm_rc_prod**3 * mm_alpha_s(3._mm_wp))) dm3a_s = dm3a_s + m3a_s_prod Index: /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_intgcm.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_intgcm.F90 (revision 3951) +++ /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_intgcm.F90 (revision 3952) @@ -8,5 +8,5 @@ ! ! The module also contains three methods: - ! - mm_initialize(dt,haze_prod_pCH4,p_prod,tx_prod,rc_prod,rm,df,rho_aer,rplanet,g0,air_rad,air_mmol,cfgpath) + ! - mm_initialize(dt,call_CH4hazeprod,p_prod,tx_prod,rc_prod,rm,df,rho_aer,rplanet,g0,air_rad,air_mmol,cfgpath) ! - read_aprm(parser,sec,pp) ! - abort_program(err) @@ -29,5 +29,5 @@ CONTAINS - SUBROUTINE mm_initialize(dt,haze_prod_pCH4,p_prod,tx_prod,rc_prod,rm,df,rho_aer,rplanet,g0,air_rad,air_mmol,clouds,cfgpath) + SUBROUTINE mm_initialize(dt,call_CH4hazeprod,p_prod,tx_prod,rc_prod,rm,df,rho_aer,rplanet,g0,air_rad,air_mmol,clouds,cfgpath) !! Initialize global parameters of the model. !! @@ -46,5 +46,5 @@ REAL(kind=mm_wp), INTENT(in) :: dt ! Enable/Disable production from CH4 photolysis. - LOGICAL :: haze_prod_pCH4 + LOGICAL :: call_CH4hazeprod ! Aerosol production pressure level (Pa). REAL(kind=mm_wp), INTENT(in) :: p_prod @@ -84,4 +84,5 @@ ! Thresholds related parameters. REAL(kind=mm_wp) :: m0as_min,rcs_min,m0af_min,rcf_min + REAL(kind=mm_wp) :: m0ccn_min,drad_min ! Debug mode control flag (may print lot of stuff if enabled). LOGICAL :: wdebug @@ -105,4 +106,6 @@ m0af_min = 1e-8_mm_wp rcf_min = 1e-9_mm_wp + m0ccn_min = 1e-8_mm_wp + drad_min = 1e-9_mm_wp wdebug = .false. @@ -132,4 +135,6 @@ err = mm_check_opt(cfg_get_value(cparser,"m0af_min",m0af_min) ,m0af_min ,1e-8_mm_wp,mm_log) err = mm_check_opt(cfg_get_value(cparser,"rcf_min",rcf_min) ,rcf_min ,rm ,mm_log) + err = mm_check_opt(cfg_get_value(cparser,"m0ccn_min",m0ccn_min) ,m0ccn_min ,1e-8_mm_wp,mm_log) + err = mm_check_opt(cfg_get_value(cparser,"drad_min",drad_min) ,drad_min ,1e-9_mm_wp,mm_log) err = mm_check_opt(cfg_get_value(cparser,"debug",wdebug) ,wdebug ,.false. ,mm_log) @@ -154,8 +159,8 @@ ! Transfert probabilities (S --> F): ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mm_w_ps2s), mm_w_ps2s, wlog=mm_log) + err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mm_call_ps2s), mm_call_ps2s, wlog=mm_log) IF (err/=0) call abort_program(err) - IF (mm_w_haze_coag .AND. mm_w_ps2s) THEN + IF (mm_call_hazecoag .AND. mm_call_ps2s) THEN err = mm_check_opt(cfg_get_value(cparser, "ps2s_file", pssfile), pssfile) IF (err /= 0) call abort_program(err) @@ -177,8 +182,8 @@ ! Mean electric correction: ! ~~~~~~~~~~~~~~~~~~~~~~~~~ - err = mm_check_opt(cfg_get_value(cparser, "electric_charging", mm_w_qe), mm_w_qe, wlog=mm_log) + err = mm_check_opt(cfg_get_value(cparser, "electric_charging", mm_call_qe), mm_call_qe, wlog=mm_log) IF (err/=0) call abort_program(err) - IF (mm_w_haze_coag .AND. mm_w_qe) THEN + IF (mm_call_hazecoag .AND. mm_call_qe) THEN err = mm_check_opt(cfg_get_value(cparser, "mq_file", mqfile), mqfile) IF (err /= 0) call abort_program(err) @@ -224,15 +229,16 @@ ! 2. YAMMS initialization !------------------------ - err = mm_global_init_0(dt,df,rm,rho_aer,haze_prod_pCH4,p_prod,tx_prod,rc_prod, & - rplanet,g0,air_rad,air_mmol,coag_choice, & - w_h_prod,w_h_sed,w_h_coag,fwsed_m0,fwsed_m3, & - m0as_min,rcs_min,m0af_min,rcf_min, & - clouds,spcpath,wdebug) + err = mm_global_init_0(dt,df,rm,rho_aer,call_CH4hazeprod,p_prod,tx_prod,rc_prod, & + rplanet,g0,air_rad,air_mmol,coag_choice, & + w_h_prod,w_h_sed,w_h_coag,fwsed_m0,fwsed_m3, & + m0as_min,rcs_min,m0af_min,rcf_min, & + clouds,spcpath,m0ccn_min,drad_min & + wdebug) IF (err /= 0) call abort_program(err) ! Dump parameters. WRITE(*,'(a)') "========= MUPHYS PARAMETERS ===========" - WRITE(*,'(a,L2)') "transfert_probability: ", mm_w_ps2s - WRITE(*,'(a,L2)') "electric_charging : ", mm_w_qe + WRITE(*,'(a,L2)') "transfert_probability: ", mm_call_ps2s + WRITE(*,'(a,L2)') "electric_charging : ", mm_call_qe call mm_dump_parameters() Index: /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_methods.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_methods.F90 (revision 3951) +++ /trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_methods.F90 (revision 3952) @@ -98,5 +98,5 @@ res = 1._mm_wp - IF (rcs <= 0.0_mm_wp .OR. .NOT.mm_w_ps2s) RETURN + IF (rcs <= 0.0_mm_wp .OR. .NOT.mm_call_ps2s) RETURN SELECT CASE(k+flow) @@ -148,5 +148,5 @@ chx = 0 - IF (.NOT.mm_w_qe) THEN + IF (.NOT.mm_call_qe) THEN res = 1._mm_wp RETURN Index: /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_diagnostics.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_diagnostics.F90 (revision 3951) +++ /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_diagnostics.F90 (revision 3952) @@ -10,7 +10,13 @@ ! ! The module also contains three methods: -! - inimuphy_diag(ngrid,nlayer,nq,pq,int2ext) -! - alpha_s(k) -! - alpha_f(k) +! - inimuphy_diag(ngrid,nlayer,nq,pq,int2ext) | Initialize the variables associated to microphysics diagnostics. +! - alpha_s(k) | Inter-moment relation for spherical aerosols size distribution law. +! - alpha_f(k) | Inter-moment relation for fractal aerosols size distribution law. +! +! @warning +! We suppose a given order of microphysical tracers in micro_indx and micro_ice_indx: +! 1. mu_m0as, 2. mu_m3as, 3. mu_m0af, 4. mu_m3af. +! If clouds are activated: +! 5. mu_m0ccn, 6. mu_m3ccn, 7(+). mu_m3ices. ! ! Authors @@ -20,7 +26,9 @@ USE tracer_h - + USE callkeys_mod, only : haze_rc_prod, haze_rm, callmuclouds + IMPLICIT NONE + ! Haze model diagnostics REAL(kind=8), ALLOCATABLE, DIMENSION(:) :: mp2m_aer_s_prec ! Spherical aerosols precipitation (kg.m-2.s-1). REAL(kind=8), ALLOCATABLE, DIMENSION(:) :: mp2m_aer_f_prec ! Fractal aerosols precipitation (kg.m-2.s-1). @@ -31,7 +39,19 @@ REAL(kind=8), ALLOCATABLE, DIMENSION(:,:) :: mp2m_rc_sph ! Spherical mode characteristic radius, i.e. bulk radius (m). REAL(kind=8), ALLOCATABLE, DIMENSION(:,:) :: mp2m_rc_fra ! Fractal mode characteristic radius, i.e. bulk radius (m). - - !$OMP THREADPRIVATE(mp2m_aer_s_prec,mp2m_aer_f_prec,mp2m_aer_s_w,mp2m_aer_f_w,mp2m_aer_s_flux,mp2m_aer_f_flux) - !$OMP THREADPRIVATE(mp2m_rc_sph,mp2m_rc_fra) +!$OMP THREADPRIVATE(mp2m_aer_s_prec,mp2m_aer_f_prec,mp2m_aer_s_w,mp2m_aer_f_w,mp2m_aer_s_flux,mp2m_aer_f_flux) +!$OMP THREADPRIVATE(mp2m_rc_sph,mp2m_rc_fra) + + ! Cloud model diagnostics + REAL(kind=8), ALLOCATABLE, DIMENSION(:) :: mp2m_ccn_prec ! CCN precipitation (kg.m-2.s-1). + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:) :: mp2m_ice_prec ! Ice precipitation (kg.m-2.s-1). + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:) :: mp2m_cld_w ! Cloud drop settling velocity (m.s-1). + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:) :: mp2m_ccn_flux ! CCN mass flux (kg.m-2.s-1). + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:,:) :: mp2m_ice_fluxes ! Ice mass fluxes (kg.m-2.s-1). + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:) :: mp2m_rc_cld ! Cloud drop radius (m). + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:,:) :: mp2m_gas_sat ! Condensible gaz saturation ratios. + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:,:) :: mp2m_nrate ! Condensible species nucleation rate (m-2.s-1). + REAL(kind=8), ALLOCATABLE, DIMENSION(:,:,:) :: mp2m_grate ! Condensible species growth rate (m2.s-1). +!$OMP THREADPRIVATE(mp2m_ccn_prec,mp2m_ice_prec,mp2m_cld_w,mp2m_ccn_flux,mp2m_ice_fluxes) +!$OMP THREADPRIVATE(mp2m_rc_cld,mp2m_gas_sat,mp2m_nrate,mp2m_grate) CONTAINS @@ -40,6 +60,4 @@ !! Initialize the variables associated to microphysics diagnostics. !! - use callkeys_mod, only : haze_rc_prod, haze_rm - INTEGER, INTENT(in) :: ngrid ! Number of points of the horizontal grid. INTEGER, INTENT(in) :: nlayer ! Number of points of the vertical grid (layers). @@ -51,8 +69,12 @@ ! Local variables: !~~~~~~~~~~~~~~~~~ - REAL, DIMENSION(:,:), ALLOCATABLE :: m0as ! 0th order moment of the spherical mode (m-3). - REAL, DIMENSION(:,:), ALLOCATABLE :: m3as ! 3rd order moment of the spherical mode (m3.m-3). - REAL, DIMENSION(:,:), ALLOCATABLE :: m0af ! 0th order moment of the fractal mode (m-3). - REAL, DIMENSION(:,:), ALLOCATABLE :: m3af ! 3rd order moment of the fractal mode (m3.m-3). + INTEGER :: i + REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m0as ! 0th order moment of the spherical mode (m-3). + REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m3as ! 3rd order moment of the spherical mode (m3.m-3). + REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m0af ! 0th order moment of the fractal mode (m-3). + REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m3af ! 3rd order moment of the fractal mode (m3.m-3). + REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m0ccn ! 0th order moment of the cloud condensation nuclei (m-3). + REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m3ccn ! 3rd order moment of the cloud condensation nuclei (m3.m-3). + REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m3itot ! 3rd order moment of the total ices (m3.m-3). IF (.NOT.ALLOCATED(mp2m_aer_s_prec)) ALLOCATE(mp2m_aer_s_prec(ngrid)) @@ -64,5 +86,25 @@ IF (.NOT.ALLOCATED(mp2m_rc_sph)) ALLOCATE(mp2m_rc_sph(ngrid,nlayer)) IF (.NOT.ALLOCATED(mp2m_rc_fra)) ALLOCATE(mp2m_rc_fra(ngrid,nlayer)) - + + IF (.NOT.ALLOCATED(mp2m_ccn_prec)) ALLOCATE(mp2m_ccn_prec(ngrid)) + IF (.NOT.ALLOCATED(mp2m_ice_prec)) ALLOCATE(mp2m_ice_prec(ngrid,nmicro_ices)) + IF (.NOT.ALLOCATED(mp2m_cld_w)) ALLOCATE(mp2m_cld_w(ngrid,nlayer)) + IF (.NOT.ALLOCATED(mp2m_ccn_flux)) ALLOCATE(mp2m_ccn_flux(ngrid,nlayer)) + IF (.NOT.ALLOCATED(mp2m_ice_fluxes)) ALLOCATE(mp2m_ice_fluxes(ngrid,nlayer,nmicro_ices)) + IF (.NOT.ALLOCATED(mp2m_rc_cld)) ALLOCATE(mp2m_rc_cld(ngrid,nlayer)) + IF (.NOT.ALLOCATED(mp2m_gas_sat)) ALLOCATE(mp2m_gas_sat(ngrid,nlayer,nmicro_ices)) + IF (.NOT.ALLOCATED(mp2m_nrate)) ALLOCATE(mp2m_nrate(ngrid,nlayer,nmicro_ices)) + IF (.NOT.ALLOCATED(mp2m_grate)) ALLOCATE(mp2m_grate(ngrid,nlayer,nmicro_ices)) + + ALLOCATE(m0as(ngrid,nlayer)) + ALLOCATE(m3as(ngrid,nlayer)) + ALLOCATE(m0af(ngrid,nlayer)) + ALLOCATE(m3af(ngrid,nlayer)) + ALLOCATE(m0ccn(ngrid,nlayer)) + ALLOCATE(m3ccn(ngrid,nlayer)) + ALLOCATE(m3itot(ngrid,nlayer)) + + ! Set to 0: + !~~~~~~~~~~ mp2m_aer_s_prec(:) = 0d0 mp2m_aer_f_prec(:) = 0d0 @@ -71,11 +113,31 @@ mp2m_aer_s_flux(:,:) = 0d0 mp2m_aer_f_flux(:,:) = 0d0 - - ! Initializes the radius to their first value (useful for radiative transfer). - m0as = pq(:,:,micro_indx(1)) * int2ext - m3as = pq(:,:,micro_indx(2)) * int2ext - m0af = pq(:,:,micro_indx(3)) * int2ext - m3af = pq(:,:,micro_indx(4)) * int2ext - + mp2m_rc_sph(:,:) = 0d0 + mp2m_rc_fra(:,:) = 0d0 + + mp2m_ccn_prec(:) = 0d0 + mp2m_ice_prec(:,:) = 0d0 + mp2m_cld_w(:,:) = 0d0 + mp2m_ccn_flux(:,:) = 0d0 + mp2m_ice_fluxes(:,:,:) = 0d0 + mp2m_rc_cld(:,:) = 0d0 + mp2m_gas_sat(:,:,:) = 0d0 + mp2m_nrate(:,:,:) = 0d0 + mp2m_grate(:,:,:) = 0d0 + + m0as(:,:) = 0d0 + m3as(:,:) = 0d0 + m0af(:,:) = 0d0 + m3af(:,:) = 0d0 + m0ccn(:,:) = 0d0 + m3ccn(:,:) = 0d0 + m3itot(:,:) = 0d0 + + ! Initializes the radius to their first value (useful for radiative transfer): + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! Spherical aerosols + m0as(:,:) = pq(:,:,micro_indx(1)) * int2ext + m3as(:,:) = pq(:,:,micro_indx(2)) * int2ext + WHERE(m0as > 1e-8 .AND. m3as > (1e-8*alpha_s(3.)*haze_rc_prod**3)) mp2m_rc_sph = (m3as / (m0as*alpha_s(3.)))**(1./3.) @@ -84,4 +146,8 @@ ENDWHERE + ! Fractal aerosols + m0af(:,:) = pq(:,:,micro_indx(3)) * int2ext + m3af(:,:) = pq(:,:,micro_indx(4)) * int2ext + WHERE(m0af > 1e-8 .AND. m3af > (1e-8*alpha_f(3.)*haze_rm**3)) mp2m_rc_fra = (m3af / (m0af*alpha_f(3.)))**(1./3.) @@ -89,5 +155,20 @@ mp2m_rc_fra = 0d0 ENDWHERE - + + ! Cloud drops + IF (callmuclouds) THEN + m0ccn(:,:) = pq(:,:,micro_indx(5)) * int2ext + m3ccn(:,:) = pq(:,:,micro_indx(6)) * int2ext + DO i=1, size(micro_ice_indx) + m3itot(:,:) = m3itot(:,:) + (pq(:,:,micro_ice_indx(i)) * int2ext) + ENDDO + + WHERE(m0ccn > 1e-8 .AND. (m3ccn+m3itot) > (1e-8*alpha_f(3.)*haze_rm**3)) + mp2m_rc_cld = ((m3ccn+m3itot) / (m0ccn))**(1./3.) + ELSEWHERE + mp2m_rc_cld = 0d0 + ENDWHERE + ENDIF ! end of callmuclouds + CONTAINS Index: /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90 (revision 3951) +++ /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90 (revision 3952) @@ -1,5 +1,5 @@ subroutine inimufi(ptimestep) use callkeys_mod, only : call_haze_prod_pCH4, haze_p_prod, haze_tx_prod, haze_rc_prod, & - haze_rm, haze_df, haze_rho, air_rad, & + haze_rm, haze_df, haze_rho, air_rad, & callmuclouds use tracer_h Index: /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3951) +++ /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3952) @@ -181,6 +181,6 @@ ! pdv(ngrid,nlayer) \ Temporal derivative of the corresponding ! pdt(ngrid,nlayer) / variables due to physical processes. -! pdq(ngrid,nlayer) / -! pdpsrf(ngrid) / +! pdq(ngrid,nlayer,nq) / +! pdpsrf(ngrid) / ! ! @@ -2528,4 +2528,19 @@ call write_output("m3"//TRIM(str(6:)),"Volume of "//TRIM(str(6:))//" ice","m3.m-3",zq(:,:,micro_ice_indx(iq))*int2ext(:,:)) enddo + + ! Diagnostics: + call write_output("ccn_prec","Cloud condensation nuclei precipitation","kg.m-2.s-1",mp2m_ccn_prec(:)) + call write_output("cld_w","Cloud drop settling velocity","m.s-1",mp2m_cld_w(:,:)) + call write_output("ccn_flux","Cloud condensation nuclei mass flux","kg.m-2.s-1",mp2m_ccn_flux(:,:)) + call write_output("rcld","Cloud drop radius","m",mp2m_rc_cld(:,:)) + + do iq = 1, size(micro_ice_indx) + str = TRIM(nameOfTracer(micro_ice_indx(iq))) + call write_output(TRIM(str(6:))//"ice_prec",TRIM(str(6:))//" ice precipitation","kg.m-2.s-1",mp2m_ice_prec(:,iq)) + call write_output(TRIM(str(6:))//"ice_flux",TRIM(str(6:))//" ice mass flux","kg.m-2.s-1",mp2m_ice_fluxes(:,:,iq)) + call write_output(TRIM(str(6:))//"_sat",TRIM(str(6:))//" saturation ratio","",mp2m_gas_sat(:,:,iq)) + call write_output(TRIM(str(6:))//"_nrate",TRIM(str(6:))//" nucleation rate","m-2.s-1",mp2m_nrate(:,:,iq)) + call write_output(TRIM(str(6:))//"_grate",TRIM(str(6:))//" ice growth rate","m2.s-1",mp2m_grate(:,:,iq)) + enddo endif ! end callmuclouds endif ! end callmufi