Changeset 3083 for trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90
- Timestamp:
- Oct 12, 2023, 10:30:22 AM (15 months ago)
- File:
-
- 1 edited
-
trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90 (modified) (44 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90
r2109 r3083 1 ! Copyright 2013-2015,2017 Université de Reims Champagne-Ardenne1 ! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne 2 2 ! Contributor: J. Burgalat (GSMA, URCA) 3 3 ! email of the author : jeremie.burgalat@univ-reims.fr 4 ! 4 ! 5 5 ! This software is a computer program whose purpose is to compute 6 6 ! microphysics processes using a two-moments scheme. 7 ! 7 ! 8 8 ! This library is governed by the CeCILL-B license under French law and 9 ! abiding by the rules of distribution of free software. You can use, 9 ! abiding by the rules of distribution of free software. You can use, 10 10 ! modify and/ or redistribute the software under the terms of the CeCILL-B 11 11 ! license as circulated by CEA, CNRS and INRIA at the following URL 12 ! "http://www.cecill.info". 13 ! 12 ! "http://www.cecill.info". 13 ! 14 14 ! As a counterpart to the access to the source code and rights to copy, 15 15 ! modify and redistribute granted by the license, users are provided only 16 16 ! with a limited warranty and the software's author, the holder of the 17 17 ! economic rights, and the successive licensors have only limited 18 ! liability. 19 ! 18 ! liability. 19 ! 20 20 ! In this respect, the user's attention is drawn to the risks associated 21 21 ! with loading, using, modifying and/or developing or reproducing the … … 25 25 ! professionals having in-depth computer knowledge. Users are therefore 26 26 ! encouraged to load and test the software's suitability as regards their 27 ! requirements in conditions enabling the security of their systems and/or 28 ! data to be ensured and, more generally, to use and operate it in the 29 ! same conditions as regards security. 30 ! 27 ! requirements in conditions enabling the security of their systems and/or 28 ! data to be ensured and, more generally, to use and operate it in the 29 ! same conditions as regards security. 30 ! 31 31 ! The fact that you are presently reading this means that you have had 32 32 ! knowledge of the CeCILL-B license and that you accept its terms. … … 35 35 !! summary: Haze microphysics module. 36 36 !! author: J. Burgalat 37 !! date: 2013-2015,2017 37 !! date: 2013-2015,2017,2022 38 38 MODULE MM_HAZE 39 39 !! Haze microphysics module. … … 46 46 !! 47 47 !! @note 48 !! The production function is specific to Titan, where aerosols are created above the detached 49 !! haze layer. No other source is taken into account. This process is controled by two parameters, 50 !! the pressure level of production and the production rate. Then both M0 and M3 of the aerosols 48 !! The production function is specific to Titan, where aerosols are created above the detached 49 !! haze layer. No other source is taken into account. This process is controled by two parameters, 50 !! the pressure level of production and the production rate. Then both M0 and M3 of the aerosols 51 51 !! distribution are updated in the production zone by addition of the production rate along a 52 52 !! gaussian shape. 53 53 !! 54 54 !! @note 55 !! The interface methods always uses the global variables defined in [[mm_globals(module)]] when 55 !! The interface methods always uses the global variables defined in [[mm_globals(module)]] when 56 56 !! values (any kind, temperature, pressure, moments...) over the vertical grid are required. 57 57 !! 58 58 !! @warning 59 !! The tendencies returned by the method are always defined over the vertical grid from __TOP__ 59 !! The tendencies returned by the method are always defined over the vertical grid from __TOP__ 60 60 !! to __GROUND__. 61 61 !! 62 !! @todo 62 !! @todo 63 63 !! Modify tests on tendencies vectors to get sure that allocation is done: 64 64 !! Currently, we assume the compiler handles automatic allocation of arrays. … … 72 72 73 73 PUBLIC :: mm_haze_microphysics, mm_haze_coagulation, mm_haze_sedimentation, & 74 75 76 74 mm_haze_production 75 76 CONTAINS 77 77 78 78 !============================================================================ … … 83 83 !! Get the evolution of moments tracers through haze microphysics processes. 84 84 !! 85 !! The subroutine is a wrapper to the haze microphysics methods. It computes the tendencies 86 !! of moments tracers for coagulation, sedimentation and production processes for the 85 !! The subroutine is a wrapper to the haze microphysics methods. It computes the tendencies 86 !! of moments tracers for coagulation, sedimentation and production processes for the 87 87 !! atmospheric column. 88 88 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_s 89 89 !! Tendency of the 0th order moment of the spherical mode distribution (\(m^{-3}\)). 90 90 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_s 91 91 !! Tendency of the 3rd order moment of the spherical mode distribution (\(m^{3}.m^{-3}\)). 92 92 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_f 93 93 !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-3}\)). 94 94 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_f 95 95 !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-3}\)). 96 96 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0as 97 97 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3as … … 107 107 zdm3af(1:mm_nla) = 0._mm_wp 108 108 109 IF (mm_w_haze_coag) THEN 109 IF (mm_w_haze_coag) THEN 110 110 ! Calls coagulation 111 111 call mm_haze_coagulation(dm0a_s,dm3a_s,dm0a_f,dm3a_f) … … 120 120 121 121 ! Updates tendencies 122 dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as 122 dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as 123 123 dm0a_f=dm0a_f+zdm0af ; dm3a_f=dm3a_f+zdm3af 124 124 ENDIF … … 127 127 call mm_haze_production(zdm0as,zdm3as) 128 128 ! We only produce spherical aerosols 129 dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as 129 dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as 130 130 ENDIF 131 131 … … 137 137 ! COAGULATION PROCESS RELATED METHODS 138 138 !============================================================================ 139 139 140 140 SUBROUTINE mm_haze_coagulation(dM0s,dM3s,dM0f,dM3f) 141 141 !! Get the evolution of the aerosols moments vertical column due to coagulation process. 142 !! 142 !! 143 143 !! This is main method of the coagulation process: 144 144 !! … … 149 149 !! 5. Finally computes the tendencies of the moments. 150 150 !! 151 !! All arguments are assumed vectors of __N__ elements where __N__ is the total number of 151 !! All arguments are assumed vectors of __N__ elements where __N__ is the total number of 152 152 !! vertical __layers__. 153 153 !! 154 154 !! @note 155 !! The method uses directly the global variables related to the vertical atmospheric structure 156 !! stored in [[mm_globals(module)]]. Consequently they must be updated before calling the subroutine. 155 !! The method uses directly the global variables related to the vertical atmospheric structure 156 !! stored in [[mm_globals(module)]]. Consequently they must be updated before calling the subroutine. 157 157 !! 158 158 !! @bug 159 !! If the transfert probabilities are set to 1 for the two flow regimes (pco and pfm), 159 !! If the transfert probabilities are set to 1 for the two flow regimes (pco and pfm), 160 160 !! a floating point exception occured (i.e. a NaN) as we perform a division by zero 161 161 !! … … 163 163 !! Get rid of the fu\*\*\*\* STOP statement... 164 164 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0s 165 165 !! Tendency of the 0th order moment of the spherical size-distribution over a time step (\(m^{-3}\)). 166 166 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3s 167 167 !! Tendency of the 3rd order moment of the spherical size-distribution (\(m^{3}.m^{-3}\)). 168 168 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0f 169 169 !! Tendency of the 0th order moment of the fractal size-distribution over a time step (\(m^{-3}\)). 170 170 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3f 171 171 !! Tendency of the 3rd order moment of the fractal size-distribution over a time step (\(m^{3}.m^{-3}\)). 172 172 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: c_kco,c_kfm,c_slf,tmp, & 173 173 kco,kfm,pco,pfm,mq 174 174 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a_ss,a_sf,b_ss,b_ff,c_ss,c_sf 175 175 INTEGER :: i … … 180 180 ! Alloctes local arrays 181 181 ALLOCATE(kco(mm_nla),kfm(mm_nla),c_slf(mm_nla), & 182 183 182 c_kco(mm_nla),c_kfm(mm_nla),mq(mm_nla), & 183 pco(mm_nla),pfm(mm_nla)) 184 184 ALLOCATE(a_ss(mm_nla),a_sf(mm_nla), & 185 186 187 185 b_ss(mm_nla),b_ff(mm_nla), & 186 c_ss(mm_nla),c_sf(mm_nla)) 187 188 188 a_ss(:) = 0._mm_wp ; a_sf(:) = 0._mm_wp 189 b_ss(:) = 0._mm_wp ; b_ff(:) = 0._mm_wp 189 b_ss(:) = 0._mm_wp ; b_ff(:) = 0._mm_wp 190 190 c_ss(:) = 0._mm_wp ; c_sf(:) = 0._mm_wp 191 191 … … 193 193 c_kco(:) = mm_get_kco(mm_temp) ; c_kfm(:) = mm_get_kfm(mm_temp) 194 194 ! get slf (slip-flow factor) 195 c_slf(:) = mm_akn * mm_lambda_g(mm_temp,mm_play) 195 c_slf(:) = mm_akn * mm_lambda_g(mm_temp,mm_play) 196 196 197 197 DO i=1,mm_nla … … 202 202 pfm(i) = mm_ps2s(mm_rcs(i),0,1,mm_temp(i),mm_play(i)) 203 203 ! (A_SS_CO x A_SS_FM) / (A_SS_CO + A_SS_FM) 204 kco(i) = g0ssco(mm_rcs(i),c_slf(i),c_kco(i)) 205 kfm(i) = g0ssfm(mm_rcs(i),c_kfm(i)) 204 kco(i) = g0ssco(mm_rcs(i),c_slf(i),c_kco(i)) 205 kfm(i) = g0ssfm(mm_rcs(i),c_kfm(i)) 206 206 IF (kco(i)*(pco(i)-2._mm_wp)+kfm(i)*(pfm(i)-2._mm_wp) /=0) THEN 207 207 a_ss(i) = (kco(i)*(pco(i)-2._mm_wp)*kfm(i)*(pfm(i)-2._mm_wp))/(kco(i)*(pco(i)-2._mm_wp)+kfm(i)*(pfm(i)-2._mm_wp)) … … 267 267 ENDIF 268 268 ENDDO 269 269 270 270 DEALLOCATE(kco,kfm,c_kco,c_kfm,pco,pfm,c_slf) 271 271 … … 302 302 !! Get γ pre-factor for the 0th order moment with SS interactions in the continuous flow regime. 303 303 !! 304 !! @note 304 !! @note 305 305 !! If __rcs__ is 0, the function returns 0. 306 306 REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. … … 320 320 !! Get γ pre-factor for the 0th order moment with SF interactions in the continuous flow regime. 321 321 !! 322 !! @note 322 !! @note 323 323 !! If __rcs__ or __rcf__ is 0, the function returns 0. 324 324 REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. … … 329 329 REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, e, rcff 330 330 res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN 331 e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra 332 ! computes mm_alpha coefficients 333 a1=mm_alpha_s(1._mm_wp) ; a2=mm_alpha_f(-e) ; a3=mm_alpha_f(e) 331 e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra 332 ! computes mm_alpha coefficients 333 a1=mm_alpha_s(1._mm_wp) ; a2=mm_alpha_f(-e) ; a3=mm_alpha_f(e) 334 334 a4=mm_alpha_s(-1._mm_wp) ; a5=mm_alpha_s(-2._mm_wp) ; a6=mm_alpha_f(-2._mm_wp*e) 335 335 ! Computes gamma pre-factor 336 336 res = c_kco*( 2._mm_wp + a1*a2*rcs/rcff + a4*a3*rcff/rcs + c_slf*( a4/rcs + & 337 337 a2/rcff + a5*a3*rcff/rcs**2 + a1*a6*rcs/rcff**2)) 338 338 RETURN 339 339 END FUNCTION g0sfco … … 342 342 !! Get γ pre-factor for the 0th order moment with FF interactions in the continuous flow regime. 343 343 !! 344 !! @note 344 !! @note 345 345 !! If __rcf__ is 0, the function returns 0. 346 346 REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. … … 361 361 !! Get γ pre-factor for the 3rd order moment with SS interactions in the continuous flow regime. 362 362 !! 363 !! @note 363 !! @note 364 364 !! If __rcs__ is 0, the function returns 0. 365 365 REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. … … 370 370 res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN 371 371 ! computes mm_alpha coefficients 372 a1=mm_alpha_s(3._mm_wp) ; a2=mm_alpha_s(2._mm_wp) ; a3=mm_alpha_s(1._mm_wp) 372 a1=mm_alpha_s(3._mm_wp) ; a2=mm_alpha_s(2._mm_wp) ; a3=mm_alpha_s(1._mm_wp) 373 373 a4=mm_alpha_s(4._mm_wp) ; a5=mm_alpha_s(-1._mm_wp) ; a6=mm_alpha_s(-2._mm_wp) 374 374 375 375 ! Computes gamma pre-factor 376 res = (2._mm_wp*a1 + a2*a3 + a4*a5 + c_slf/rcs*(a3**2 + a4*a6 + a1*a5 + a2))* & 377 376 res = (2._mm_wp*a1 + a2*a3 + a4*a5 + c_slf/rcs*(a3**2 + a4*a6 + a1*a5 + a2))* & 377 c_kco/(a1**2*rcs**3) 378 378 RETURN 379 379 END FUNCTION g3ssco … … 382 382 !! Get γ pre-factor for the 3rd order moment with SF interactions in the continuous flow regime. 383 383 !! 384 !! @note 384 !! @note 385 385 !! If __rcs__ or __rcf__ is 0, the function returns 0. 386 386 REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. … … 392 392 res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN 393 393 ! computes mm_alpha coefficients 394 e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra 395 a1=mm_alpha_s(3._mm_wp) ; a2=mm_alpha_s(4._mm_wp) ; a3=mm_alpha_f(-e) 394 e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra 395 a1=mm_alpha_s(3._mm_wp) ; a2=mm_alpha_s(4._mm_wp) ; a3=mm_alpha_f(-e) 396 396 a4=mm_alpha_s(2._mm_wp) ; a5=mm_alpha_f(e) ; a6=mm_alpha_s(1._mm_wp) 397 397 a7=mm_alpha_f(-2._mm_wp*e) ; a8=mm_alpha_f(3._mm_wp) 398 398 ! Computes gamma pre-factor 399 399 res = (2._mm_wp*a1*rcs**3 + a2*rcs**4*a3/rcff + a4*rcs**2*a5*rcff + & 400 401 400 c_slf *( a4*rcs**2 + a1*rcs**3*a3/rcff + a6*rcs*a5*rcff + & 401 a2*rcs**4*a7/rcff**2))* c_kco/(a1*a8*(rcs*rcf)**3) 402 402 RETURN 403 403 END FUNCTION g3sfco … … 406 406 !! Get γ pre-factor for the 0th order moment with SS interactions in the Free Molecular flow regime. 407 407 !! 408 !! @note 408 !! @note 409 409 !! If __rcs__ is 0, the function returns 0. 410 410 REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. … … 422 422 423 423 ELEMENTAL FUNCTION g0sffm(rcs, rcf, c_kfm) RESULT(res) 424 !> Get γ pre-factor for the 0th order moment with SF interactions in the Free Molecular flow regime. 425 !! 426 !! @note 424 !> Get γ pre-factor for the 0th order moment with SF interactions in the Free Molecular flow regime. 425 !! 426 !! @note 427 427 !! If __rcs__ or __rcf__ is 0, the function returns 0. 428 428 REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. … … 435 435 res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN 436 436 ! computes mm_alpha coefficients 437 e1 = 3._mm_wp/mm_df 438 e2 = 6._mm_wp/mm_df 437 e1 = 3._mm_wp/mm_df 438 e2 = 6._mm_wp/mm_df 439 439 e3 = (6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 440 440 e4 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 441 441 442 rcff1 = mm_rb2ra * rcf**e1 442 rcff1 = mm_rb2ra * rcf**e1 443 443 rcff2 = rcff1**2 444 rcff3 = mm_rb2ra * rcf**e3 444 rcff3 = mm_rb2ra * rcf**e3 445 445 rcff4 = mm_rb2ra**2 * rcf**e4 446 446 447 a1=mm_alpha_s(0.5_mm_wp) ; a2=mm_alpha_s(-0.5_mm_wp) ; a3=mm_alpha_f(e1) 447 a1=mm_alpha_s(0.5_mm_wp) ; a2=mm_alpha_s(-0.5_mm_wp) ; a3=mm_alpha_f(e1) 448 448 a4=mm_alpha_s(-1.5_mm_wp) ; a5=mm_alpha_f(e2) ; a6=mm_alpha_s(2._mm_wp) 449 a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(1._mm_wp) ; a9=mm_alpha_f(e3) 449 a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(1._mm_wp) ; a9=mm_alpha_f(e3) 450 450 a10=mm_alpha_f(e4) 451 451 452 452 ! Computes gamma pre-factor 453 453 res = (a1*rcs**0.5_mm_wp + 2._mm_wp*rcff1*a2*a3/rcs**0.5_mm_wp + a4*a5*rcff2/rcs**1.5_mm_wp + & 454 455 454 a6*a7*rcs**2/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs*rcff3 + a10*rcff4 & 455 )*mm_get_btk(4,0)*c_kfm 456 456 RETURN 457 457 END FUNCTION g0sffm 458 458 459 459 ELEMENTAL FUNCTION g0fffm(rcf, c_kfm) RESULT(res) 460 !! Get γ pre-factor for the 0th order moment with FF interactions in the Free Molecular flow regime. 461 !! 462 !! @note 460 !! Get γ pre-factor for the 0th order moment with FF interactions in the Free Molecular flow regime. 461 !! 462 !! @note 463 463 !! If __rcf__ is 0, the function returns 0. 464 464 REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. 465 465 REAL(kind=mm_wp), INTENT(in) :: c_kfm !! Thermodynamic free molecular flow regime pre-factor. 466 REAL(kind=mm_wp) :: res !! γ coagulation kernel pre-factor. 466 REAL(kind=mm_wp) :: res !! γ coagulation kernel pre-factor. 467 467 REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, e1, e2, e3, rcff 468 468 res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN 469 469 ! computes mm_alpha coefficients 470 e1=3._mm_wp/mm_df ; e2=(6_mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 470 e1=3._mm_wp/mm_df ; e2=(6_mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 471 471 e3=(12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 472 472 rcff=mm_rb2ra**2*rcf**e3 … … 489 489 res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN 490 490 ! computes mm_alpha coefficients 491 a1=mm_alpha_s(3.5_mm_wp) ; a2=mm_alpha_s(1._mm_wp) ; a3=mm_alpha_s(2.5_mm_wp) 491 a1=mm_alpha_s(3.5_mm_wp) ; a2=mm_alpha_s(1._mm_wp) ; a3=mm_alpha_s(2.5_mm_wp) 492 492 a4=mm_alpha_s(2._mm_wp) ; a5=mm_alpha_s(1.5_mm_wp) ; a6=mm_alpha_s(5._mm_wp) 493 a7=mm_alpha_s(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp) ; a9=mm_alpha_s(-0.5_mm_wp) 493 a7=mm_alpha_s(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp) ; a9=mm_alpha_s(-0.5_mm_wp) 494 494 a10=mm_alpha_s(3._mm_wp) ; a11=mm_alpha_s(0.5_mm_wp) 495 495 ! Computes gamma pre-factor 496 496 res = (a1 + 2._mm_wp*a2*a3 + a4*a5 + a6*a7 + 2._mm_wp*a8*a9 + a10*a11) & 497 497 *mm_get_btk(1,3)*c_kfm/(a10**2*rcs**2.5_mm_wp) 498 498 RETURN 499 499 END FUNCTION g3ssfm … … 502 502 !! Get γ pre-factor for the 3rd order moment with SF interactions in the Free Molecular flow regime. 503 503 !! 504 !! @note 504 !! @note 505 505 !! If __rcs__ or __rcf__ is 0, the function returns 0. 506 506 REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. … … 512 512 res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN 513 513 ! computes mm_alpha coefficients 514 e1=3._mm_wp/mm_df 515 e2=(6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 514 e1=3._mm_wp/mm_df 515 e2=(6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 516 516 e3=(12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 517 517 rcff1=mm_rb2ra*rcf**e1 ; rcff2=mm_rb2ra*rcf**e2 ; rcff3=mm_rb2ra**2*rcf**e3 518 a1=mm_alpha_s(3.5_mm_wp) ; a2=mm_alpha_s(2.5_mm_wp) ; a3=mm_alpha_f(e1) 519 a4=mm_alpha_s(1.5_mm_wp) ; a5=mm_alpha_f(2._mm_wp*e1) ; a6=mm_alpha_s(5._mm_wp) 520 a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp) ; a9=mm_alpha_f(e2) 518 a1=mm_alpha_s(3.5_mm_wp) ; a2=mm_alpha_s(2.5_mm_wp) ; a3=mm_alpha_f(e1) 519 a4=mm_alpha_s(1.5_mm_wp) ; a5=mm_alpha_f(2._mm_wp*e1) ; a6=mm_alpha_s(5._mm_wp) 520 a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp) ; a9=mm_alpha_f(e2) 521 521 a10=mm_alpha_s(3._mm_wp) ; a11=mm_alpha_f(e3) ; a12=mm_alpha_f(3._mm_wp) 522 522 ! Computes gamma pre-factor 523 523 res = (a1*rcs**3.5_mm_wp + 2._mm_wp*a2*a3*rcs**2.5_mm_wp*rcff1 + a4*a5*rcs**1.5_mm_wp*rcff1**2 + & 524 a6*a7*rcs**5/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs**4*rcff2 + &525 524 a6*a7*rcs**5/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs**4*rcff2 + & 525 a10*a11*rcs**3*rcff3)*mm_get_btk(4,3)*c_kfm/(a10*a12*(rcs*rcf)**3) 526 526 RETURN 527 527 END FUNCTION g3sffm … … 530 530 ! SEDIMENTATION PROCESS RELATED METHODS 531 531 !============================================================================ 532 532 533 533 SUBROUTINE mm_haze_sedimentation(dm0s,dm3s,dm0f,dm3f) 534 534 !! Interface to sedimentation algorithm. 535 535 !! 536 536 !! The subroutine computes the evolution of each moment of the aerosols tracers 537 !! through sedimentation process and returns their tendencies for a timestep. 537 !! through sedimentation process and returns their tendencies for a timestep. 538 538 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s 539 539 !! Tendency of the 0th order moment of the spherical mode (\(m^{-3}\)). 540 540 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s 541 541 !! Tendency of the 3rd order moment of the spherical mode (\(m^{3}.m^{-3}\)). 542 542 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0f 543 543 !! Tendency of the 0th order moment of the fractal mode (\(m^{-3}\)). 544 544 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3f 545 545 !! Tendency of the 3rd order moment of the fractal mode (\(m^{3}.m^{-3}\)). 546 546 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: ft,fdcor,wth 547 REAL(kind=mm_wp) :: m,n,p 548 REAL(kind=mm_wp), PARAMETER :: fac = 4._mm_wp/3._mm_wp * mm_pi 547 REAL(kind=mm_wp), PARAMETER :: fac = 4._mm_wp/3._mm_wp * mm_pi 549 548 550 549 ALLOCATE(ft(mm_nle),wth(mm_nle),fdcor(mm_nle)) … … 624 623 !! Compute the tendency of the moment through sedimentation process. 625 624 !! 626 !! 625 !! 627 626 !! The method computes the time evolution of the \(k^{th}\) order moment through sedimentation: 628 627 !! … … 630 629 !! 631 630 !! The equation is resolved using a [Crank-Nicolson algorithm](http://en.wikipedia.org/wiki/Crank-Nicolson_method). 632 !! 633 !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells 631 !! 632 !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells 634 633 !! from ancient times. It is based on \cite{toon1988b,fiadeiro1977,turco1979} and is an update of the algorithm 635 634 !! originally implemented in the LMDZ-Titan 2D GCM. … … 638 637 REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: ft !! Downward sedimentation flux (effective velocity of the moment). 639 638 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dmk !! Tendency of \(k^{th}\) order moment (in \(m^{k}.m^{-3}\)). 640 INTEGER :: i 639 INTEGER :: i 641 640 REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: as,bs,cs,mko 642 641 ALLOCATE(as(mm_nla), bs(mm_nla), cs(mm_nla), mko(mm_nla)) … … 648 647 bs(1:mm_nla) = -(ft(2:mm_nle)+mm_dzlay(1:mm_nla)/dt) 649 648 cs(1:mm_nla) = -mm_dzlay(1:mm_nla)/dt*mk(1:mm_nla) 650 ! (Tri)diagonal matrix inversion 649 ! (Tri)diagonal matrix inversion 651 650 mko(1) = cs(1)/bs(1) 652 651 DO i=2,mm_nla ; mko(i) = (cs(i)-mko(i-1)*as(i))/bs(i) ; ENDDO … … 660 659 ! interior 661 660 mko(2:mm_nla-1)=(bs(2:mm_nla-1)*mk(1:mm_nla-2) + & 662 663 661 cs(2:mm_nla-1)*mk(2:mm_nla-1) & 662 )/as(2:mm_nla-1) 664 663 ENDIF 665 664 dmk = mko - mk … … 670 669 SUBROUTINE get_weff(mk,k,df,rc,dt,afun,wth,corf) 671 670 !! Get the effective settling velocity for aerosols moments. 672 !! 673 !! This method computes the effective settling velocity of the \(k^{th}\) order moment of aerosol 674 !! tracers. The basic settling velocity (\(v^{eff}_{M_{k}}\)) is computed using the following 671 !! 672 !! This method computes the effective settling velocity of the \(k^{th}\) order moment of aerosol 673 !! tracers. The basic settling velocity (\(v^{eff}_{M_{k}}\)) is computed using the following 675 674 !! equation: 676 !! 677 !! $$ 675 !! 676 !! $$ 678 677 !! \begin{eqnarray*} 679 !! \Phi^{sed}_{M_{k}} &=& \int_{0}^{\infty} n(r) r^{k} \times w(r) dr 678 !! \Phi^{sed}_{M_{k}} &=& \int_{0}^{\infty} n(r) r^{k} \times w(r) dr 680 679 !! == M_{k} \times v^{eff}_{M_{k}} \\ 681 680 !! v^{eff}_{M_{k} &=& \dfrac{2 \rho g r_{c}^{\dfrac{3D_{f}-3}{D_{f}}}} … … 686 685 !! $$ 687 686 !! 688 !! \(v^{eff}_{M_{k}\) is then corrected to reduce numerical diffusion of the sedimentation algorithm 687 !! \(v^{eff}_{M_{k}\) is then corrected to reduce numerical diffusion of the sedimentation algorithm 689 688 !! as defined in \cite{toon1988b}. 690 689 !! 691 690 !! @warning 692 !! Both __df__, __rc__ and __afun__ must be consistent with each other otherwise wrong values will 691 !! Both __df__, __rc__ and __afun__ must be consistent with each other otherwise wrong values will 693 692 !! be computed. 694 693 REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla) :: mk 695 694 !! Moment of order __k__ (\(m^{k}.m^{-3}\)) at each layer. 696 695 REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla) :: rc 697 696 !! Characteristic radius associated to the moment at each layer. 698 697 REAL(kind=mm_wp), INTENT(in) :: k 699 698 !! The order of the moment. 700 699 REAL(kind=mm_wp), INTENT(in) :: df 701 700 !! Fractal dimension of the aersols. 702 701 REAL(kind=mm_wp), INTENT(in) :: dt 703 702 !! Time step (s). 704 703 REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle) :: wth 705 706 REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle), OPTIONAL :: corf 707 704 !! Theoretical Settling velocity at each vertical __levels__ (\( wth \times corf = weff\)). 705 REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle), OPTIONAL :: corf 706 !! _Fiadero_ correction factor applied to the theoretical settling velocity at each vertical __levels__. 708 707 INTERFACE 709 708 FUNCTION afun(order) … … 713 712 REAL(kind=mm_wp) :: afun !! Alpha value. 714 713 END FUNCTION afun 715 END INTERFACE 714 END INTERFACE 716 715 INTEGER :: i 717 716 REAL(kind=mm_wp) :: af1,af2,ar1,ar2 … … 720 719 REAL(kind=mm_wp), DIMENSION(mm_nle) :: zcorf 721 720 ! ------------------ 722 721 723 722 wth(:) = 0._mm_wp ; zcorf(:) = 1._mm_wp 724 723 725 724 ar1 = (3._mm_wp*df -3._mm_wp)/df ; ar2 = (3._mm_wp*df -6._mm_wp)/df 726 af1 = (df*(k+3._mm_wp)-3._mm_wp)/df ; af2 = (df*(k+3._mm_wp)-6._mm_wp)/df 725 af1 = (df*(k+3._mm_wp)-3._mm_wp)/df ; af2 = (df*(k+3._mm_wp)-6._mm_wp)/df 727 726 rb2ra = mm_rm**((df-3._mm_wp)/df) 728 727 DO i=2,mm_nla 729 IF (rc(i-1) <= 0._mm_wp) CYCLE 728 IF (rc(i-1) <= 0._mm_wp) CYCLE 730 729 dzb = (mm_dzlay(i)+mm_dzlay(i-1))/2._mm_wp 731 730 csto = 2._mm_wp*mm_rhoaer*mm_effg(mm_zlev(i))/(9._mm_wp*mm_eta_g(mm_btemp(i))) 732 731 cslf = mm_akn * mm_lambda_g(mm_btemp(i),mm_plev(i)) 733 wth(i) = - csto/(rb2ra*afun(k)) * (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2)) 732 wth(i) = - csto/(rb2ra*afun(k)) * (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2)) 733 734 ! >>> [TEMPO : BBT] 735 !wth(i) = wth(i) * (2574e3 / (2574e3+mm_zlev(i)))**4 736 ! <<< [TEMPO : BBT] 737 734 738 ! now correct velocity to reduce numerical diffusion 735 739 IF (.NOT.mm_no_fiadero_w) THEN 736 740 IF (mk(i) <= 0._mm_wp) THEN 737 ratio = mm_fiadero_max 741 ratio = mm_fiadero_max 738 742 ELSE 739 ratio = MAX(MIN(mk(i-1)/mk(i),mm_fiadero_max),mm_fiadero_min) 743 ratio = MAX(MIN(mk(i-1)/mk(i),mm_fiadero_max),mm_fiadero_min) 740 744 ENDIF 741 745 ! apply correction … … 760 764 SUBROUTINE mm_haze_production(dm0s,dm3s) 761 765 !! Compute the production of aerosols moments. 762 !! 763 !! The method computes the tendencies of M0 and M3 for the spherical mode through production process. 764 !! Production values are distributed along a normal law in altitude, centered at 766 !! 767 !! The method computes the tendencies of M0 and M3 for the spherical mode through production process. 768 !! Production values are distributed along a normal law in altitude, centered at 765 769 !! [[mm_globals(module):mm_p_prod(variable)]] pressure level with a fixed sigma of 20km. 766 770 !! 767 !! First M3 tendency is computed and M0 is retrieved using the inter-moments relation a spherical 771 !! First M3 tendency is computed and M0 is retrieved using the inter-moments relation a spherical 768 772 !! characteristic radius set to [[mm_globals(module):mm_rc_prod(variable)]]. 769 773 !! … … 773 777 REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s !! Tendency of M3 (\(m^{3}.m^{-3}\)). 774 778 INTEGER :: i 775 REAL(kind=mm_wp) :: zprod,cprod,timefact 779 REAL(kind=mm_wp) :: zprod,cprod,timefact 776 780 REAL(kind=mm_wp), PARAMETER :: sigz = 20e3_mm_wp, & 777 781 fnorm = 1._mm_wp/(dsqrt(2._mm_wp*mm_pi)*sigz), & … … 793 797 dm3s(:)= mm_tx_prod *0.75_mm_wp/mm_pi *mm_dt / mm_rhoaer / 2._mm_wp / mm_dzlev(1:mm_nla) * & 794 798 (erf((mm_zlev(1:mm_nla)-zprod)/znorm) - & 795 erf((mm_zlev(2:mm_nla+1)-zprod)/znorm)) 799 erf((mm_zlev(2:mm_nla+1)-zprod)/znorm)) 796 800 dm0s(:) = dm3s(:)/(mm_rc_prod**3*mm_alpha_s(3._mm_wp)) 797 801 … … 803 807 ENDIF 804 808 805 806 809 END SUBROUTINE mm_haze_production 807 810
Note: See TracChangeset
for help on using the changeset viewer.
