From 307ecdfd01a8e17e0d228de1833d0d667bb52c22 Mon Sep 17 00:00:00 2001 From: Philippe Le Sager Date: Fri, 29 May 2026 15:24:46 +0200 Subject: [PATCH 1/2] Add remark from Xuemei regarding test on volume in M&N --- ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 b/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 index 4bdee7c..2824879 100644 --- a/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 +++ b/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 @@ -527,6 +527,12 @@ SUBROUTINE AER_ACTIV_MORALES_NENES_FULL(KIDIA, KFDIA, KTDIA, KLON, KLEV, LCLOUD, NCL(JL) = NNACL(JL) NH2SO4(JL) = NSO4(JL) - NNA2SO4(JL) + ! QUESTION: should we use a fixed lower value (e.g. + ! 1e-30) to avoid ovestimating aerosol volume by SP. + ! + ! The argument is that we miss all value between + ! EPS(SP)=~1e-7 and EPS(DP)=~2e-16 when running at SP + ! IF (ZVOL(JL) .GE. ZEPS) THEN !eehol: total volume per mode need to be above treshold to avoid div by zero !---mode kappa = volume-weighted sum of component kappa's ZKAPPA(JL,JK,JMOD) = ( (Kap_ss * NNACL(JL) * WNACL / (DNACL*1.E3_JPRB)) + & From afe4d713c47961f5dc5b55fd847fed52d0f839e4 Mon Sep 17 00:00:00 2001 From: Philippe Le Sager Date: Fri, 29 May 2026 16:06:49 +0200 Subject: [PATCH 2/2] Use the standard erf implementation to compute CDF --- ifs-source/arpifs/m7/module/mo_ham_m7.F90 | 198 ++-------------------- 1 file changed, 11 insertions(+), 187 deletions(-) diff --git a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 index ad94206..afee439 100644 --- a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 @@ -59,16 +59,15 @@ MODULE mo_ham_m7 -USE mo_kind, ONLY: dp + USE mo_kind, ONLY: dp -IMPLICIT NONE - -PRIVATE + IMPLICIT NONE -PUBLIC :: m7, m7_cumulative_normal + PRIVATE + PUBLIC :: m7, m7_cumulative_normal -!REAL(dp), PUBLIC, ALLOCATABLE :: rwet_m7(:,:,:), rdry_m7(:,:,:) -!REAL(dp), PUBLIC, ALLOCATABLE :: densaer_m7(:,:,:), aerwat_m7(:,:,:) + !REAL(dp), PUBLIC, ALLOCATABLE :: rwet_m7(:,:,:), rdry_m7(:,:,:) + !REAL(dp), PUBLIC, ALLOCATABLE :: densaer_m7(:,:,:), aerwat_m7(:,:,:) CONTAINS @@ -143,192 +142,17 @@ SUBROUTINE m7_cumulative_normal ( arg, presult, ccum ) ! functions, and proper selection of the machine-dependent ! constants. ! - ! Explanation of machine-dependent constants. - ! - ! MIN = smallest machine representable number. - ! - ! EPS = argument below which anorm(x) may be represented by - ! 0.5 and above which x*x will not underflow. - ! A conservative value is the largest machine number X - ! such that 1.0 + X = 1.0 to machine precision. - ! - ! Error returns - ! - ! The program returns ANORM = 0 for ARG .LE. XLOW. - ! - ! Author: - ! - ! W. J. Cody - ! Mathematics and Computer Science Division - ! Argonne National Laboratory - ! Argonne, IL 60439 - ! - ! Latest modification: March 15, 1992 ! USE mo_kind, ONLY: dp ! IMPLICIT NONE ! - REAL(dp), PARAMETER, DIMENSION ( 5 ) :: a = (/ & - 2.2352520354606839287e00_dp, & - 1.6102823106855587881e02_dp, & - 1.0676894854603709582e03_dp, & - 1.8154981253343561249e04_dp, & - 6.5682337918207449113e-2_dp /) - REAL(dp) :: arg - REAL(dp), PARAMETER, DIMENSION ( 4 ) :: b = (/ & - 4.7202581904688241870e01_dp, & - 9.7609855173777669322e02_dp, & - 1.0260932208618978205e04_dp, & - 4.5507789335026729956e04_dp /) - REAL(dp), PARAMETER, DIMENSION ( 9 ) :: c = (/ & - 3.9894151208813466764e-1_dp, & - 8.8831497943883759412e00_dp, & - 9.3506656132177855979e01_dp, & - 5.9727027639480026226e02_dp, & - 2.4945375852903726711e03_dp, & - 6.8481904505362823326e03_dp, & - 1.1602651437647350124e04_dp, & - 9.8427148383839780218e03_dp, & - 1.0765576773720192317e-8_dp /) - REAL(dp) :: ccum - REAL(dp), PARAMETER, DIMENSION ( 8 ) :: d = (/ & - 2.2266688044328115691e01_dp, & - 2.3538790178262499861e02_dp, & - 1.5193775994075548050e03_dp, & - 6.4855582982667607550e03_dp, & - 1.8615571640885098091e04_dp, & - 3.4900952721145977266e04_dp, & - 3.8912003286093271411e04_dp, & - 1.9685429676859990727e04_dp /) - REAL(dp) :: del -!@@@ REAL(dp) :: dpmpar - REAL(dp) :: eps - INTEGER :: i - REAL(dp) :: zmin - REAL(dp), PARAMETER, DIMENSION ( 6 ) :: p = (/ & - 2.1589853405795699e-1_dp, & - 1.274011611602473639e-1_dp, & - 2.2235277870649807e-2_dp, & - 1.421619193227893466e-3_dp, & - 2.9112874951168792e-5_dp, & - 2.307344176494017303e-2_dp /) - REAL(dp), PARAMETER, DIMENSION ( 5 ) :: q = (/ & - 1.28426009614491121e00_dp, & - 4.68238212480865118e-1_dp, & - 6.59881378689285515e-2_dp, & - 3.78239633202758244e-3_dp, & - 7.29751555083966205e-5_dp /) - REAL(dp) :: presult - REAL(dp), PARAMETER :: root32 = 5.656854248_dp - REAL(dp), PARAMETER :: sixten = 16.0_dp - REAL(dp) :: temp - REAL(dp), PARAMETER :: sqrpi = 3.9894228040143267794e-1_dp - REAL(dp), PARAMETER :: thrsh = 0.66291_dp - REAL(dp) :: x - REAL(dp) :: xden - REAL(dp) :: xnum - REAL(dp) :: y - REAL(dp) :: xsq - - - !REAL(dp), ALLOCATABLE :: rwet_m7(:,:,:), rdry_m7(:,:,:) - !REAL(dp), ALLOCATABLE :: densaer_m7(:,:,:), aerwat_m7(:,:,:) - ! - ! Machine dependent constants - ! - eps = EPSILON ( 1.0_dp ) * 0.5_dp - ! - !@@@ Simplified calculation of the smallest machine representable number - ! (Higher accuracy than needed!) - ! - !@@@ min = dpmpar(2) - - zmin = EPSILON ( 1.0_dp ) - - x = arg - y = ABS ( x ) - - IF ( y <= thrsh ) THEN - ! - ! Evaluate anorm for |X| <= 0.66291 - ! - IF ( y > eps ) THEN - xsq = x * x - ELSE - xsq = 0.0_dp - END IF + REAL(dp), INTENT(IN) :: arg + REAL(dp), INTENT(OUT) :: ccum + REAL(dp), INTENT(OUT) :: presult - xnum = a(5) * xsq - xden = xsq - DO i = 1, 3 - xnum = ( xnum + a(i) ) * xsq - xden = ( xden + b(i) ) * xsq - END DO - presult = x * ( xnum + a(4) ) / ( xden + b(4) ) - temp = presult - presult = 0.5_dp + temp - ccum = 0.5_dp - temp - ! - ! Evaluate ANORM for 0.66291 <= |X| <= sqrt(32) - ! - ELSE IF ( y <= root32 ) THEN - - xnum = c(9) * y - xden = y -!DIR$ UNROLL -!CDIR UNROLL=7 - DO i = 1, 7 - xnum = ( xnum + c(i) ) * y - xden = ( xden + d(i) ) * y - END DO - presult = ( xnum + c(8) ) / ( xden + d(8) ) - xsq = AINT ( y * sixten ) / sixten - del = ( y - xsq ) * ( y + xsq ) - presult = EXP(-xsq*xsq*0.5_dp) * EXP(-del*0.5_dp) * presult - ccum = 1.0_dp - presult - - IF ( x > 0.0_dp ) THEN - temp = presult - presult = ccum - ccum = temp - END IF - ! - ! Evaluate anorm for |X| > sqrt(32). - ! - ELSE - - presult = 0.0_dp - xsq = 1.0_dp / ( x * x ) - xnum = p(6) * xsq - xden = xsq - DO i = 1, 4 - xnum = ( xnum + p(i) ) * xsq - xden = ( xden + q(i) ) * xsq - END DO - - presult = xsq * ( xnum + p(5) ) / ( xden + q(5) ) - presult = ( sqrpi - presult ) / y - xsq = AINT ( x * sixten ) / sixten - del = ( x - xsq ) * ( x + xsq ) - presult = EXP ( - xsq * xsq * 0.5_dp ) * EXP ( - del * 0.5_dp ) * presult - ccum = 1.0_dp - presult - - IF ( x > 0.0_dp ) THEN - temp = presult - presult = ccum - ccum = temp - END IF - - END IF - - IF ( presult < zmin ) THEN - presult = 0.0_dp - END IF - - IF ( ccum < zmin ) THEN - ccum = 0.0_dp - END IF + presult = 0.5_dp * (1.0_dp + erf(arg / sqrt(2.0_dp))) + ccum = 1.0_dp - presult END SUBROUTINE m7_cumulative_normal