Changeset 4032 for trunk/LMDZ.PLUTO/libf/muphypluto
- Timestamp:
- Jan 28, 2026, 11:32:11 AM (4 weeks ago)
- Location:
- trunk/LMDZ.PLUTO/libf/muphypluto
- Files:
-
- 2 edited
-
mp2m_globals.F90 (modified) (4 diffs)
-
mp2m_haze.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_globals.F90
r3957 r4032 30 30 ! - Vertical structure part 31 31 ! 32 ! The module also contains twenty methods:32 ! The module also contains twenty-one methods: 33 33 ! - mm_global_init_0 34 34 ! - mm_global_init_1 … … 36 36 ! - mm_aerosols_init 37 37 ! - mm_clouds_init 38 ! - mm_alpha_s, mm_alpha_f 38 ! - mm_alpha_s, mm_alpha_f, mm_alpha_ccn 39 39 ! - mm_set_moments_thresholds 40 40 ! - mm_get_rcs, mm_get_rcf … … 80 80 ! Moments parameters (mm_aerosols_init) 81 81 PROTECTED :: mm_m0aer_s, mm_m3aer_s, mm_m0aer_f, mm_m3aer_f 82 ! Moments parameters (derived, are updated with moments parameters)83 PROTECTED :: mm_rcs, mm_rcf84 82 ! Thresholds parameters 85 83 PROTECTED :: mm_m0as_min, mm_m3as_min, mm_rcs_min, mm_m0af_min, mm_m3af_min, mm_rcf_min … … 1170 1168 END FUNCTION mm_alpha_f 1171 1169 1170 PURE FUNCTION mm_alpha_ccn(k) RESULT (res) 1171 !! Inter-moment relation for fractal aerosols size distribution law. 1172 !! Mk / M0 = rc^k . alpha(k) 1173 !! 1174 REAL(kind=mm_wp), INTENT(in) :: k ! k order of the moment. 1175 REAL(kind=mm_wp) :: sigma ! Standard deviation. 1176 REAL(kind=mm_wp) :: res ! Alpha value. 1177 1178 ! Pluto's case 1179 !~~~~~~~~~~~~~ 1180 sigma = 0.2_mm_wp 1181 res = exp(k**2 * sigma**2 / 2._mm_wp) 1182 1183 RETURN 1184 END FUNCTION mm_alpha_ccn 1185 1172 1186 1173 1187 !============================================================================ -
trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_haze.F90
r3952 r4032 76 76 ! Local variables: 77 77 !~~~~~~~~~~~~~~~~~ 78 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0as 79 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3as 80 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0af 81 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3af 78 ! Updated tracers with haze related processes tendencies 79 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: m0as,m3as 80 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: m0af,m3af 81 ! Temporary tendencies through production/coagulation/sedimentation 82 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: pdm0as,pdm3as 83 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: cdm0as,cdm3as,cdm0af,cdm3af 84 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0as,zdm3as,zdm0af,zdm3af 82 85 83 86 ! Initialization: 84 87 !~~~~~~~~~~~~~~~~ 85 dm0a_s = 0._mm_wp 86 dm3a_s = 0._mm_wp 87 dm0a_f = 0._mm_wp 88 dm3a_f = 0._mm_wp 89 88 ALLOCATE(m0as(mm_nla),m3as(mm_nla),m0af(mm_nla),m3af(mm_nla)) 89 ALLOCATE(pdm0as(mm_nla),pdm3as(mm_nla)) 90 ALLOCATE(cdm0as(mm_nla),cdm3as(mm_nla),cdm0af(mm_nla),cdm3af(mm_nla)) 90 91 ALLOCATE(zdm0as(mm_nla),zdm3as(mm_nla),zdm0af(mm_nla),zdm3af(mm_nla)) 91 zdm0as(1:mm_nla) = 0._mm_wp 92 zdm3as(1:mm_nla) = 0._mm_wp 93 zdm0af(1:mm_nla) = 0._mm_wp 94 zdm3af(1:mm_nla) = 0._mm_wp 95 92 93 m0as(:) = mm_m0aer_s ; m3as(:) = mm_m3aer_s 94 m0af(:) = mm_m0aer_f ; m3af(:) = mm_m3aer_f 95 96 dm0a_s = 0._mm_wp ; dm3a_s = 0._mm_wp 97 dm0a_f = 0._mm_wp ; dm3a_f = 0._mm_wp 98 99 pdm0as(1:mm_nla) = 0._mm_wp ; pdm3as(1:mm_nla) = 0._mm_wp 100 cdm0as(1:mm_nla) = 0._mm_wp ; cdm3as(1:mm_nla) = 0._mm_wp 101 cdm0af(1:mm_nla) = 0._mm_wp ; cdm3af(1:mm_nla) = 0._mm_wp 102 zdm0as(1:mm_nla) = 0._mm_wp ; zdm3as(1:mm_nla) = 0._mm_wp 103 zdm0af(1:mm_nla) = 0._mm_wp ; zdm3af(1:mm_nla) = 0._mm_wp 96 104 97 105 ! Microphysical processes: 98 106 !~~~~~~~~~~~~~~~~~~~~~~~~~ 99 107 100 ! Coagulation101 IF (mm_call_hazecoag) THEN102 call mm_haze_coagulation(dm0a_s,dm3a_s,dm0a_f,dm3a_f)103 ENDIF104 105 ! Sedimentation106 IF (mm_call_hazesed) THEN107 call mm_haze_sedimentation(zdm0as,zdm3as,zdm0af,zdm3af)108 109 ! Computes precipitation110 mm_aers_prec = SUM(zdm3as*mm_dzlev*mm_rhoaer)111 mm_aerf_prec = SUM(zdm3af*mm_dzlev*mm_rhoaer)112 113 ! Updates tendencies114 dm0a_s = dm0a_s + zdm0as115 dm3a_s = dm3a_s + zdm3as116 dm0a_f = dm0a_f + zdm0af117 dm3a_f = dm3a_f + zdm3af118 ENDIF119 120 108 ! Production 121 109 IF (mm_call_hazeprod) THEN 122 110 ! Production by photolysis of CH4 123 111 IF (mm_call_CH4hazeprod) THEN 124 dm0a_s = dm0a_s + (m3a_s_prod / (mm_rc_prod**3 * mm_alpha_s(3._mm_wp)))125 dm3a_s = dm3a_s +m3a_s_prod112 pdm0as = m3a_s_prod / (mm_rc_prod**3 * mm_alpha_s(3._mm_wp)) 113 pdm3as = m3a_s_prod 126 114 ELSE 127 call mm_haze_production(zdm0as,zdm3as) 128 dm0a_s = dm0a_s + zdm0as 129 dm3a_s = dm3a_s + zdm3as 115 call mm_haze_production(pdm0as,pdm3as) 130 116 ENDIF 117 118 ! Update tracers 119 m0as = m0as + pdm0as ; m3as = m3as + pdm3as 120 WHERE(m3as > 0._mm_wp .AND. m0as > 0._mm_wp) 121 mm_rcs = mm_get_rcs(m0as,m3as) 122 ELSEWHERE 123 mm_rcs = 0._mm_wp 124 ENDWHERE 131 125 ENDIF 126 127 ! Coagulation 128 IF (mm_call_hazecoag) THEN 129 call mm_haze_coagulation(m0as,m3as,m0af,m3af,cdm0as,cdm3as,cdm0af,cdm3af) 130 131 ! Update tracers 132 m0as = m0as + cdm0as ; m3as = m3as + cdm3as 133 m0af = m0af + cdm0af ; m3af = m3af + cdm3af 134 WHERE(m3as > 0._mm_wp .AND. m0as > 0._mm_wp) 135 mm_rcs = mm_get_rcs(m0as,m3as) 136 ELSEWHERE 137 mm_rcs = 0._mm_wp 138 ENDWHERE 139 WHERE(m3af > 0._mm_wp .AND. m0af > 0._mm_wp) 140 mm_rcf = mm_get_rcf(m0af,m3af) 141 ELSEWHERE 142 mm_rcf = 0._mm_wp 143 ENDWHERE 144 ENDIF 145 146 ! Sedimentation 147 IF (mm_call_hazesed) THEN 148 call mm_haze_sedimentation(m0as,m3as,m0af,m3af,zdm0as,zdm3as,zdm0af,zdm3af) 149 150 ! Update tracers 151 m0as = m0as + zdm0as ; m3as = m3as + zdm3as 152 m0af = m0af + zdm0af ; m3af = m3af + zdm3af 153 WHERE(m3as > 0._mm_wp .AND. m0as > 0._mm_wp) 154 mm_rcs = mm_get_rcs(m0as,m3as) 155 ELSEWHERE 156 mm_rcs = 0._mm_wp 157 ENDWHERE 158 WHERE(m3af > 0._mm_wp .AND. m0af > 0._mm_wp) 159 mm_rcf = mm_get_rcf(m0af,m3af) 160 ELSEWHERE 161 mm_rcf = 0._mm_wp 162 ENDWHERE 163 164 ! Computes precipitation 165 mm_aers_prec = SUM(zdm3as*mm_dzlev*mm_rhoaer) 166 mm_aerf_prec = SUM(zdm3af*mm_dzlev*mm_rhoaer) 167 ENDIF 168 169 ! Updates tendencies: 170 !~~~~~~~~~~~~~~~~~~~~ 171 dm0a_s = pdm0as + cdm0as + zdm0as 172 dm3a_s = pdm3as + cdm3as + zdm3as 173 dm0a_f = cdm0af + zdm0af 174 dm3a_f = cdm3af + zdm3af 132 175 133 176 RETURN … … 187 230 !============================================================================ 188 231 189 SUBROUTINE mm_haze_sedimentation( dm0s,dm3s,dm0f,dm3f)232 SUBROUTINE mm_haze_sedimentation(m0as,m3as,m0af,m3af,dm0s,dm3s,dm0f,dm3f) 190 233 !! Interface to sedimentation algorithm. 191 234 !! … … 193 236 !! through sedimentation process and returns their tendencies for a timestep. 194 237 !! 238 ! 0th and 3rd order moment of the aerosols (spherical mode) component. 239 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0as ! (m-3) 240 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3as ! (m3.m-3) 241 ! 0th and 3rd order moment of the aerosols (fractal mode) component. 242 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0af ! (m-3) 243 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3af ! (m3.m-3) 195 244 ! Tendency of the 0th order moment of the spherical mode (m-3). 196 245 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s … … 214 263 call get_weff(0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) 215 264 216 call let_me_fall_in_peace(m m_m0aer_s,-wth,mm_dt,dm0s)217 call let_me_fall_in_peace(m m_m3aer_s,-wth,mm_dt,dm3s)265 call let_me_fall_in_peace(m0as,-wth,mm_dt,dm0s) 266 call let_me_fall_in_peace(m3as,-wth,mm_dt,dm3s) 218 267 219 268 ! Get settling velocity and mass flux 220 269 mm_m0as_vsed(:) = wth(1:mm_nla) 221 270 mm_m3as_vsed(:) = wth(1:mm_nla) 222 mm_aer_s_flux(:) = (m m_m3aer_s * mm_rhoaer) * wth(2:)271 mm_aer_s_flux(:) = (m3as * mm_rhoaer) * wth(2:) 223 272 224 273 ! Fractal particles 225 274 call get_weff(0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) 226 275 227 call let_me_fall_in_peace(m m_m0aer_f,-wth,mm_dt,dm0f)228 call let_me_fall_in_peace(m m_m3aer_f,-wth,mm_dt,dm3f)276 call let_me_fall_in_peace(m0af,-wth,mm_dt,dm0f) 277 call let_me_fall_in_peace(m3af,-wth,mm_dt,dm3f) 229 278 230 279 ! Get settling velocity and mass flux 231 280 mm_m0af_vsed(:) = wth(1:mm_nla) 232 281 mm_m3af_vsed(:) = wth(1:mm_nla) 233 mm_aer_f_flux(:) = (m m_m3aer_f * mm_rhoaer) * wth(2:)282 mm_aer_f_flux(:) = (m3af * mm_rhoaer) * wth(2:) 234 283 235 284 ! Force settling velocity to M3 … … 239 288 call get_weff(3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) 240 289 241 call let_me_fall_in_peace(m m_m3aer_s,-wth,mm_dt,dm3s)242 call let_me_fall_in_peace(m m_m0aer_s,-wth,mm_dt,dm0s)290 call let_me_fall_in_peace(m3as,-wth,mm_dt,dm3s) 291 call let_me_fall_in_peace(m0as,-wth,mm_dt,dm0s) 243 292 244 293 ! Get settling velocity and mass flux 245 294 mm_m3as_vsed(:) = wth(1:mm_nla) 246 295 mm_m0as_vsed(:) = wth(1:mm_nla) 247 mm_aer_s_flux(:) = (m m_m3aer_s * mm_rhoaer) * wth(2:)296 mm_aer_s_flux(:) = (m3as * mm_rhoaer) * wth(2:) 248 297 249 298 ! Fractal particles 250 299 call get_weff(3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) 251 300 252 call let_me_fall_in_peace(m m_m3aer_f,-wth,mm_dt,dm3f)253 call let_me_fall_in_peace(m m_m0aer_f,-wth,mm_dt,dm0f)301 call let_me_fall_in_peace(m3af,-wth,mm_dt,dm3f) 302 call let_me_fall_in_peace(m0af,-wth,mm_dt,dm0f) 254 303 255 304 ! Get settling velocity and mass flux 256 305 mm_m3af_vsed(:) = wth(1:mm_nla) 257 306 mm_m0af_vsed(:) = wth(1:mm_nla) 258 mm_aer_f_flux(:) = (m m_m3aer_f * mm_rhoaer) * wth(2:)307 mm_aer_f_flux(:) = (m3af * mm_rhoaer) * wth(2:) 259 308 260 309 ! No forcing of settling velocity (might be a problem...) … … 263 312 ! Spherical particles 264 313 call get_weff(0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) 265 call let_me_fall_in_peace(m m_m0aer_s,-wth,mm_dt,dm0s)314 call let_me_fall_in_peace(m0as,-wth,mm_dt,dm0s) 266 315 267 316 ! Get settling velocity … … 269 318 270 319 call get_weff(3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) 271 call let_me_fall_in_peace(m m_m3aer_s,-wth,mm_dt,dm3s)320 call let_me_fall_in_peace(m3as,-wth,mm_dt,dm3s) 272 321 273 322 ! Get settling velocity … … 275 324 276 325 ! Get mass flux 277 mm_aer_s_flux(:) = (m m_m3aer_s * mm_rhoaer) * wth(2:)326 mm_aer_s_flux(:) = (m3as * mm_rhoaer) * wth(2:) 278 327 279 328 ! Fractal aerosols 280 329 call get_weff(0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) 281 call let_me_fall_in_peace(m m_m0aer_f,-wth,mm_dt,dm0f)330 call let_me_fall_in_peace(m0af,-wth,mm_dt,dm0f) 282 331 283 332 ! Get settling velocity … … 285 334 286 335 call get_weff(3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) 287 call let_me_fall_in_peace(m m_m3aer_f,-wth,mm_dt,dm3f)336 call let_me_fall_in_peace(m3af,-wth,mm_dt,dm3f) 288 337 289 338 ! Get settling velocity … … 291 340 292 341 ! Get mass flux 293 mm_aer_f_flux(:) = (m m_m3aer_f * mm_rhoaer) * wth(2:)342 mm_aer_f_flux(:) = (m3af * mm_rhoaer) * wth(2:) 294 343 ENDIF 295 344 … … 457 506 !============================================================================ 458 507 459 SUBROUTINE mm_haze_coagulation( dM0s,dM3s,dM0f,dM3f)508 SUBROUTINE mm_haze_coagulation(m0as,m3as,m0af,m3af,dM0s,dM3s,dM0f,dM3f) 460 509 !! Get the evolution of the aerosols moments vertical column due to coagulation process. 461 510 !! … … 471 520 !! a floating point exception occured (i.e. a NaN) as we perform a division by zero 472 521 !! 522 ! 0th and 3rd order moment of the aerosols (spherical mode) component. 523 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0as ! (m-3) 524 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3as ! (m3.m-3) 525 ! 0th and 3rd order moment of the aerosols (fractal mode) component. 526 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0af ! (m-3) 527 REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3af ! (m3.m-3) 473 528 ! Tendency of the 0th order moment of the spherical size-distribution over a time step (m-3). 474 529 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0s … … 661 716 ! When X appears as square we set one X at t+1, the other a t. 662 717 ! --- dM0s 663 tmp(:) = mm_dt * (a_ss*m m_m0aer_s - a_sf*mm_m0aer_f)664 dm0s(:) = m m_m0aer_s * (tmp / (1._mm_wp - tmp))718 tmp(:) = mm_dt * (a_ss*m0as - a_sf*m0af) 719 dm0s(:) = m0as * (tmp / (1._mm_wp - tmp)) 665 720 ! --- dM0f 666 tmp(:) = mm_dt * b_ff * m m_m0aer_f667 dm0f(:) = (b_ss*mm_dt*m m_m0aer_s**2 - tmp*mm_m0aer_f) / (1._mm_wp + tmp)721 tmp(:) = mm_dt * b_ff * m0af 722 dm0f(:) = (b_ss*mm_dt*m0as**2 - tmp*m0af) / (1._mm_wp + tmp) 668 723 ! --- dM3s 669 tmp(:) = mm_dt * (c_ss*m m_m3aer_s - c_sf*mm_m3aer_f)670 dm3s(:) = m m_m3aer_s * (tmp / (1._mm_wp - tmp))724 tmp(:) = mm_dt * (c_ss*m3as - c_sf*m3af) 725 dm3s(:) = m3as * (tmp / (1._mm_wp - tmp)) 671 726 ! --- dM3f 672 727 dm3f(:) = -dm3s
Note: See TracChangeset
for help on using the changeset viewer.
