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/6] 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/6] 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 From 2ae9b4aaf9811fddb8b45830ade8be0cb4acd40a Mon Sep 17 00:00:00 2001 From: Philippe Le Sager Date: Wed, 3 Jun 2026 11:37:08 +0200 Subject: [PATCH 3/6] Use the same EPS in value in SP and DP across the M7 code --- ifs-source/arpifs/m7/module/mo_ham_activ.F90 | 10 ++--- ifs-source/arpifs/m7/module/mo_ham_m7.F90 | 12 +++--- ifs-source/arpifs/m7/module/mo_ham_rad.F90 | 19 +++++----- ifs-source/arpifs/m7/module/mo_ham_tools.F90 | 4 +- ifs-source/arpifs/m7/module/mo_ham_wetdep.F90 | 4 +- ifs-source/arpifs/m7/module/mo_kind.F90 | 3 ++ .../arpifs/m7/module/mo_math_constants.F90 | 4 +- .../arpifs/m7/module/mo_tracer_processes.F90 | 6 +-- ifs-source/arpifs/m7/module/nd_param.F90 | 6 ++- .../arpifs/m7/phys_ec/yoe_aer_activ.F90 | 37 ++++++------------- 10 files changed, 49 insertions(+), 56 deletions(-) diff --git a/ifs-source/arpifs/m7/module/mo_ham_activ.F90 b/ifs-source/arpifs/m7/module/mo_ham_activ.F90 index 73adbbe..ced5558 100644 --- a/ifs-source/arpifs/m7/module/mo_ham_activ.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_activ.F90 @@ -45,7 +45,7 @@ MODULE mo_ham_activ - USE mo_kind, ONLY: dp + USE mo_kind, ONLY: dp, THRESHOLD USE mo_ham, ONLY: sizeclass,nclass !#ifdef _OPENMP @@ -207,7 +207,7 @@ SUBROUTINE ham_activ_abdulrazzak_ghan(kproma, kbdim, klev, krow, ktdia, & pfracn(1:kproma,:,:) = 0._dp pcdncact(1:kproma,:) = 0._dp - zeps=EPSILON(1._dp) + zeps=THRESHOLD ! EPSILON(1._dp) !--- Conversions to SI units [g mol-1 to kg mol-1]: @@ -558,7 +558,7 @@ SUBROUTINE ham_activ_koehler_ab(kproma, kbdim, klev, krow, ktdia, & zsumtop(1:kproma,:,:) = 0._dp zsumbot(1:kproma,:,:) = 0._dp - zeps=EPSILON(1._dp) + zeps=THRESHOLD !EPSILON(1._dp) !--- Conversions to SI units [g mol-1 to kg mol-1]: @@ -696,7 +696,7 @@ SUBROUTINE ham_avail_activ_lin_leaitch(kproma, kbdim, klev, krow, & DO jclass=1, nclass !SF the ham_m7_logtail calculation is now done mode per mode for better efficiency IF (sizeclass(jclass)%lactivation) THEN - IF (cfracn(jclass) > EPSILON(1._dp)) THEN !SF #279: only performs this calculation when relevant + IF (cfracn(jclass) > THRESHOLD ) THEN !SF #279: only performs this calculation when relevant !SF stratiform: zr(1:kproma,:,jclass)=crcut @@ -786,7 +786,7 @@ SUBROUTINE ham_activ_diag_lin_leaitch(kproma, kbdim, klev, krow, prho, pxtm1, pc REAL(dp) :: zeps REAL(dp) :: ztmp1(kbdim,klev), ztmp2(kbdim,klev) - zeps = EPSILON(1._dp) + zeps = THRESHOLD !EPSILON(1._dp) DO jclass=1, nclass !>>dod #377 diff --git a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 index afee439..7e2555a 100644 --- a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 @@ -59,7 +59,7 @@ MODULE mo_ham_m7 - USE mo_kind, ONLY: dp + USE mo_kind, ONLY: dp, THRESHOLD IMPLICIT NONE @@ -2211,7 +2211,7 @@ SUBROUTINE m7_nuck(kproma, kbdim, klev, krow, & zqtmst = 1.0_dp/ztmst - zeps = EPSILON(1.0_dp) + zeps = THRESHOLD !EPSILON(1.0_dp) ! Relative humidity [%]: @@ -2863,7 +2863,7 @@ SUBROUTINE m7_dconc(kproma, kbdim, klev, krow, paerml, paernl, pm6dry) zfconm(:,:,:) = 1._dp zfconn(:,:,:) = 1._dp - zeps = EPSILON(1._dp) + zeps = THRESHOLD !EPSILON(1._dp) ! !--- 1) Identify how much the mode jclass has grown into the next higher mode @@ -3315,7 +3315,7 @@ SUBROUTINE m7_coaset(kproma, kbdim, klev, krow, paernl, ptp1, & zhu2, zhu !---executable procedure - zeps = EPSILON(1._dp) + zeps = THRESHOLD !EPSILON(1._dp) !---1) gridpoint properties DO jk=1,klev @@ -3556,7 +3556,7 @@ SUBROUTINE m7_concoag (kproma, kbdim, klev, krow, & !--- 0) Initializations: - zeps=EPSILON(1._dp) + zeps=THRESHOLD !EPSILON(1._dp) !--- 1) Redistribution of mass and numbers after nucleation, coagulation ---- @@ -3923,7 +3923,7 @@ SUBROUTINE m7_delcoa(kproma, kbdim, klev, krow, paerml, & !--- 0) Initialisations: ------------------------------------------------ ztmst = time_step_len - zeps = EPSILON(1._dp) + zeps = THRESHOLD !EPSILON(1._dp) za4av = 0._dp za4av1(:,:) = 0._dp diff --git a/ifs-source/arpifs/m7/module/mo_ham_rad.F90 b/ifs-source/arpifs/m7/module/mo_ham_rad.F90 index b2c34bf..fb2def3 100755 --- a/ifs-source/arpifs/m7/module/mo_ham_rad.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_rad.F90 @@ -68,7 +68,7 @@ MODULE mo_ham_rad sizeclass, & nclass USE mo_ham, ONLY: subm_aerospec - USE mo_kind, ONLY: dp + USE mo_kind, ONLY: dp, THRESHOLD USE mo_species, ONLY: speclist, naerospec, nmaxspec !>>dod soa USE mo_ham_species, ONLY: id_oc, id_wat !!mgs!! , naerospec=>ham_naerospec, aerospec=>ham_aerospec @@ -240,7 +240,7 @@ SUBROUTINE ham_rad_refrac_volume(kproma, kbdim, klev, krow, ktrac, kmod, kwv, & !---executable procedure - zeps=EPSILON(1.0_dp) + zeps=THRESHOLD !EPSILON(1.0_dp) !>>dod openmp bugfix removed allocation of arrays @@ -393,7 +393,7 @@ SUBROUTINE ham_rad_refrac_maxgar(kproma, kbdim, klev, krow, ktrac, kmod, kwv, & !---executable procedure - zeps=EPSILON(1.0_dp) + zeps=THRESHOLD !EPSILON(1.0_dp) !>>dod openmp bugfix removed allocation of arrays !<>dod deleted allocation of arrays !<EPSILON(1.0_dp)) THEN + ! PLS - TODO: SAFE DIVISION + IF (aer_piz_sw_vr(jl,jk,jwv) > THRESHOLD ) THEN aer_cg_sw_vr(jl,jk,jwv) =aer_cg_sw_vr(jl,jk,jwv)/aer_piz_sw_vr(jl,jk,jwv) END IF - IF(aer_tau_sw_vr(jl,jk,jwv)>EPSILON(1.0_dp)) THEN + IF (aer_tau_sw_vr(jl,jk,jwv) > THRESHOLD ) THEN aer_piz_sw_vr(jl,jk,jwv)=aer_piz_sw_vr(jl,jk,jwv)/aer_tau_sw_vr(jl,jk,jwv) END IF END DO @@ -1498,7 +1499,7 @@ SUBROUTINE ham_rad(kproma, kbdim, klev, krow, kpband, kb_sw, & DO jk=1, klev DO jl=1, kproma - IF(zaer_ssa_diag(jl,jk,jwv)>EPSILON(1.0_dp)) THEN + IF (zaer_ssa_diag(jl,jk,jwv) > THRESHOLD) THEN zaer_asym_diag(jl,jk,jwv) = zaer_asym_diag(jl,jk,jwv)/zaer_ssa_diag(jl,jk,jwv) zaer_ssa_diag(jl,jk,jwv) = zaer_ssa_diag(jl,jk,jwv)/ zaer_tau_diag(jl,jk,jwv) END IF @@ -2048,7 +2049,7 @@ SUBROUTINE ham_rad_diag(kproma, kbdim, klev, krow, pxtm1) !--- 0) - zeps=EPSILON(1.0_dp) + zeps = THRESHOLD ! EPSILON(1.0_dp) !--- Optical thickness for optional wavelengths: diff --git a/ifs-source/arpifs/m7/module/mo_ham_tools.F90 b/ifs-source/arpifs/m7/module/mo_ham_tools.F90 index 086c004..cfe4b33 100644 --- a/ifs-source/arpifs/m7/module/mo_ham_tools.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_tools.F90 @@ -41,7 +41,7 @@ MODULE mo_ham_tools ! *mo_ham_tools* hold auxiliary routines for the ! HAM aerosol model - USE mo_kind, ONLY: dp + USE mo_kind, ONLY: dp, THRESHOLD USE mo_exception, ONLY: finish IMPLICIT NONE @@ -275,7 +275,7 @@ SUBROUTINE ham_m7_logtail(kproma, kbdim, klev, krow, kmod, & !--- 0) - zeps=EPSILON(1.0_dp) + zeps= THRESHOLD !EPSILON(1.0_dp) !--- 1) diff --git a/ifs-source/arpifs/m7/module/mo_ham_wetdep.F90 b/ifs-source/arpifs/m7/module/mo_ham_wetdep.F90 index e7d4978..d2170c5 100755 --- a/ifs-source/arpifs/m7/module/mo_ham_wetdep.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_wetdep.F90 @@ -48,7 +48,7 @@ MODULE mo_ham_wetdep - USE mo_kind, ONLY: dp + USE mo_kind, ONLY: dp, THRESHOLD USE mo_physical_constants, ONLY: tmelt USE mo_exception, ONLY: finish USE mo_tracdef, ONLY: ntrac, trlist, AEROSOLNUMBER, AEROSOLMASS @@ -74,7 +74,7 @@ MODULE mo_ham_wetdep !--- Constants: REAL(dp), PARAMETER :: zmin = 1.e-10_dp REAL(dp), PARAMETER :: UNDEF = -999._dp - REAL(dp), PARAMETER :: zeps = EPSILON(1._dp) + REAL(dp), PARAMETER :: zeps = THRESHOLD !EPSILON(1._dp) REAL(dp), PARAMETER :: zeps_mass = 1.e-30_dp !--- Mode-wise scavenging coefficients and related diff --git a/ifs-source/arpifs/m7/module/mo_kind.F90 b/ifs-source/arpifs/m7/module/mo_kind.F90 index 1dc979b..10615e2 100644 --- a/ifs-source/arpifs/m7/module/mo_kind.F90 +++ b/ifs-source/arpifs/m7/module/mo_kind.F90 @@ -7,4 +7,7 @@ MODULE mo_kind USE parkind1, ONLY: JPRD IMPLICIT NONE ! <-- thk + + REAL(KIND=dp), PARAMETER :: THRESHOLD = REAL(EPSILON(1.0_JPRD), KIND=dp) + END MODULE mo_kind diff --git a/ifs-source/arpifs/m7/module/mo_math_constants.F90 b/ifs-source/arpifs/m7/module/mo_math_constants.F90 index 80a56ce..bd5a42b 100644 --- a/ifs-source/arpifs/m7/module/mo_math_constants.F90 +++ b/ifs-source/arpifs/m7/module/mo_math_constants.F90 @@ -19,7 +19,7 @@ MODULE mo_math_constants -USE mo_kind, ONLY: wp +USE mo_kind, ONLY: wp, THRESHOLD IMPLICIT NONE @@ -80,7 +80,7 @@ MODULE mo_math_constants !$ACC DECLARE COPYIN(rad2deg) REAL (wp), PARAMETER :: deg2rad = pi/180.0_wp REAL (wp), PARAMETER :: eps = 1.e-8_wp -REAL (wp), PARAMETER :: dbl_eps = EPSILON(1._wp) +REAL (wp), PARAMETER :: dbl_eps = THRESHOLD !EPSILON(1._wp) ! DOES NOT SEEM TO BE USED ANYWHERE REAL (wp), PARAMETER :: pi_180 = pi/180._wp ! diff --git a/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 b/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 index 8f6e03a..6aaf5e1 100755 --- a/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 +++ b/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 @@ -4,7 +4,7 @@ MODULE mo_tracer_processes - USE mo_kind, ONLY: wp + USE mo_kind, ONLY: wp, THRESHOLD IMPLICIT NONE @@ -101,7 +101,7 @@ SUBROUTINE xt_conv_massfix (kproma, kbdim, klev, & !--- 2) Mass fix mode: - zeps=EPSILON(1.0_wp) + zeps=THRESHOLD !EPSILON(1.0_wp) !--- 2.1) Calculate auxiliary variable dp/g : !--- Uppermost level: @@ -209,7 +209,7 @@ SUBROUTINE xt_borrow(kproma, kbdim, klev, klevp1, ktrac, & !--- 0) Initializations: ztmst=time_step_len - zeps=10._wp*EPSILON(1.0_wp) + zeps=10._wp*THRESHOLD !EPSILON(1.0_wp) !--- 1) Calculate auxiliary variable dp/g : diff --git a/ifs-source/arpifs/m7/module/nd_param.F90 b/ifs-source/arpifs/m7/module/nd_param.F90 index 5e08284..9af97b4 100644 --- a/ifs-source/arpifs/m7/module/nd_param.F90 +++ b/ifs-source/arpifs/m7/module/nd_param.F90 @@ -9,7 +9,8 @@ MODULE ND_PARAM !---Inherited functions, types, variables and constants USE PARKIND1, ONLY: JPIM, JPRB - + USE MO_KIND, ONLY: THRESHOLD + IMPLICIT NONE PRIVATE @@ -515,7 +516,8 @@ SUBROUTINE SINTEGRAL (SPAR, SUMMA, WPARCEL, XFHH, BET2, & ! ** Population Splitting -- Modified by Ricardo Morales 2014 DESCR = 1._JPRB - (16._JPRB/9._JPRB)*ALFA*WPARCEL*BET2*(AKOH/SPAR**2)**2 - IF (DESCR.LT.EPSILON(0.0_JPRB)) THEN + + IF (DESCR.LT.THRESHOLD) THEN CRIT2 = .TRUE. scrit = ((16._JPRB/9._JPRB)*ALFA*WPARCEL*BET2*(AKOH**2))**(0.25_JPRB) ! Scrit - (only for DELTA < 0 ) RATIO = (2.0E+7_JPRB/3.0_JPRB)*AKOH*(SPAR**(-0.3824_JPRB)-scrit**(-0.3824_JPRB)) ! Computing sp1 and sp2 (sp1 = sp2) 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 2824879..fa227db 100644 --- a/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 +++ b/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 @@ -377,24 +377,17 @@ SUBROUTINE AER_ACTIV_MORALES_NENES_FULL(KIDIA, KFDIA, KTDIA, KLON, KLEV, LCLOUD, ! Seinfeld and Pandis, Atmospheric Chemistry and Physics, Second Edition (referred to as SP) ! Morales and Nenes, JGR, D18220, 2010 - USE YOMCST, ONLY: RG, RPI - USE YOMLUN, ONLY: NULOUT - USE TM5M7_DATA, ONLY: NMOD, NSOL, DDUST, DNACL, & - & DOC, DBC, DH2SO4, DNA2SO4, DNH4NO3, DMSA, & - & NH4NO3_FACTOR, Kap_su,Kap_pom,Kap_soa, & - & Kap_bc,Kap_ss,Kap_du,Kap_na2so4,Kap_msa, & - & Kap_no3, WSO4, WH2SO4, WNACL, WNA2SO4, & - & WH2O, WDAIR - USE MO_HAM_M7CTL, ONLY: CMR2RAM, SIGMA, SIGMALN - ! USE YOE_AERO_M7_DATA, ONLY: NMOD, NSOL, SIGMA, SIGMALN, CMR2RAM, & - ! & DH2SO4, DBC, DOC, DNACL, DDUST, & - ! & DNA2SO4, DNH4NO3, DMSA, NH4NO3_FACTOR, & - ! & PPKAPPA_H2SO4, PPKAPPA_NACL, PPKAPPA_NA2SO4, & - ! & PPKAPPA_BC, PPKAPPA_OC, PPKAPPA_DU, & - ! & PPKAPPA_NH4NO3, PPKAPPA_MSA, & - ! & WSO4, WH2SO4, WNACL, WNA2SO4, & - ! & WH2O, WDAIR - USE ND_PARAM, ONLY: CCNSPEC, PDFACTIV, NDPARAM + USE YOMCST, ONLY: RG, RPI + USE YOMLUN, ONLY: NULOUT + USE TM5M7_DATA, ONLY: NMOD, NSOL, DDUST, DNACL, & + & DOC, DBC, DH2SO4, DNA2SO4, DNH4NO3, DMSA, & + & NH4NO3_FACTOR, Kap_su,Kap_pom,Kap_soa, & + & Kap_bc,Kap_ss,Kap_du,Kap_na2so4,Kap_msa, & + & Kap_no3, WSO4, WH2SO4, WNACL, WNA2SO4, & + & WH2O, WDAIR + USE MO_HAM_M7CTL, ONLY: CMR2RAM, SIGMA, SIGMALN + USE MO_KIND, ONLY: THRESHOLD + USE ND_PARAM, ONLY: CCNSPEC, PDFACTIV, NDPARAM IMPLICIT NONE @@ -488,7 +481,7 @@ SUBROUTINE AER_ACTIV_MORALES_NENES_FULL(KIDIA, KFDIA, KTDIA, KLON, KLEV, LCLOUD, PCDNC(KIDIA:KFDIA,KTDIA:KLEV) = 0._JPRB PSMAX(KIDIA:KFDIA,KTDIA:KLEV) = 0._JPRB - ZEPS=EPSILON(1._JPRB) + ZEPS=THRESHOLD !EPSILON(1._JPRB) ZWLARGE(KIDIA:KFDIA,KTDIA:KLEV) = -1._JPRB* PVERVEL(KIDIA:KFDIA,KTDIA:KLEV) / & & (RG*PRHO(KIDIA:KFDIA,KTDIA:KLEV)) @@ -527,12 +520,6 @@ 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 baad38c37a844a40b3180c8f68a750ab9360ca44 Mon Sep 17 00:00:00 2001 From: Philippe Le Sager Date: Wed, 3 Jun 2026 11:44:53 +0200 Subject: [PATCH 4/6] Flag use of EPS for safe division --- ifs-source/arpifs/m7/module/mo_ham_activ.F90 | 2 ++ ifs-source/arpifs/m7/module/mo_ham_m7.F90 | 1 + ifs-source/arpifs/m7/module/mo_ham_rad.F90 | 13 +++++++++---- ifs-source/arpifs/m7/module/mo_tracer_processes.F90 | 3 ++- ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 | 1 + 5 files changed, 15 insertions(+), 5 deletions(-) diff --git a/ifs-source/arpifs/m7/module/mo_ham_activ.F90 b/ifs-source/arpifs/m7/module/mo_ham_activ.F90 index ced5558..3afc98d 100644 --- a/ifs-source/arpifs/m7/module/mo_ham_activ.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_activ.F90 @@ -320,6 +320,7 @@ SUBROUTINE ham_activ_abdulrazzak_ghan(kproma, kbdim, klev, krow, ktdia, & !< zeps) psmax(jl,jk,:)=1._dp/SQRT(zsum(:)) ELSEWHERE @@ -596,6 +597,7 @@ SUBROUTINE ham_activ_koehler_ab(kproma, kbdim, klev, krow, ktdia, & jclass = aerocomp(jn)%iclass IF (nion > 0 .AND. sizeclass(jclass)%lactivation) THEN !>>dod<< #377 + ! PLS - TODO: SAFE DIVISION WHERE(zmasssum(1:kproma,:,jclass)>zeps) zmassfrac(1:kproma,:)=pxtm1(1:kproma,:,jt)/zmasssum(1:kproma,:,jclass) diff --git a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 index 7e2555a..7087165 100644 --- a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 @@ -4304,6 +4304,7 @@ SUBROUTINE m7_delcoa(kproma, kbdim, klev, krow, paerml, & DO kmod=1,nclass IF (kmod>jclass) THEN !kai + ! PLS - TODO: SAFE DIVISION IF (abs(zbftot).gt.zeps) then pbfract1(jl,jk,kmod-jclass)=pbfract1(jl,jk,kmod-jclass)/zbftot ELSE diff --git a/ifs-source/arpifs/m7/module/mo_ham_rad.F90 b/ifs-source/arpifs/m7/module/mo_ham_rad.F90 index fb2def3..26a5e5c 100755 --- a/ifs-source/arpifs/m7/module/mo_ham_rad.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_rad.F90 @@ -304,6 +304,7 @@ SUBROUTINE ham_rad_refrac_volume(kproma, kbdim, klev, krow, ktrac, kmod, kwv, & !---Weighted averaging: DO jk=1,klev DO jl=1,kproma + ! PLS - TODO: SAFE DIVISION IF(zvsum(jl,jk)>zeps) THEN pnr(jl,jk)=znrsum(jl,jk)/zvsum(jl,jk) @@ -466,7 +467,7 @@ SUBROUTINE ham_rad_refrac_maxgar(kproma, kbdim, klev, krow, ktrac, kmod, kwv, & DO jk=1,klev DO jl=1,kproma !<zeps) THEN cn_eff(jl,jk)=CMPLX( znrsum(jl,jk)/zvsum(jl,jk) , znisum(jl,jk)/zvsum(jl,jk), kind=dp ) ELSE @@ -486,7 +487,7 @@ SUBROUTINE ham_rad_refrac_maxgar(kproma, kbdim, klev, krow, ktrac, kmod, kwv, & DO jk=1,klev DO jl=1,kproma !<zeps) THEN IF (zvcore(jl,jk)/zvsum(jl,jk)>zeps .AND. zvcore(jl,jk)/zvsum(jl,jk)<(1.0_dp-zeps)) THEN lcore(jl,jk)=.TRUE. @@ -556,7 +557,7 @@ SUBROUTINE ham_rad_refrac_maxgar(kproma, kbdim, klev, krow, ktrac, kmod, kwv, & DO jk=1,klev DO jl=1,kproma - + ! PLS - TODO: SAFE DIVISION IF(zvsum(jl,jk)>zeps .AND. lcore(jl,jk)) THEN cn_0(jl,jk)=CMPLX( znrsum(jl,jk)/zvsum(jl,jk) , znisum(jl,jk)/zvsum(jl,jk), kind=dp ) ELSE @@ -1499,6 +1500,7 @@ SUBROUTINE ham_rad(kproma, kbdim, klev, krow, kpband, kb_sw, & DO jk=1, klev DO jl=1, kproma + ! PLS - TODO: SAFE DIVISION IF (zaer_ssa_diag(jl,jk,jwv) > THRESHOLD) THEN zaer_asym_diag(jl,jk,jwv) = zaer_asym_diag(jl,jk,jwv)/zaer_ssa_diag(jl,jk,jwv) zaer_ssa_diag(jl,jk,jwv) = zaer_ssa_diag(jl,jk,jwv)/ zaer_tau_diag(jl,jk,jwv) @@ -2121,6 +2123,7 @@ SUBROUTINE ham_rad_diag(kproma, kbdim, klev, krow, pxtm1) END DO !>>SF #458 (replacing WHERE statements) + ! PLS - TODO: SAFE DIVISION ll1(1:kproma) = (ztau(1:kproma) > zeps) ztmp1(1:kproma) = MERGE(ztau(1:kproma), 1._dp, ll1(1:kproma)) !SF 1._dp is a dummy val. @@ -2273,7 +2276,8 @@ SUBROUTINE ham_rad_diag(kproma, kbdim, klev, krow, pxtm1) DO jk=1, klev DO jl=1, kproma - IF (zvsum(jl,jk,jclass)>zeps) THEN + ! PLS - TODO: SAFE DIVISION + IF (zvsum(jl,jk,jclass)>zeps) THEN ztaucomp(jl)=ztaucomp(jl) + & tau_p(jl,jk,krow)*zvcomp(jl,jk,jspec,jclass)/zvsum(jl,jk,jclass) zabscomp(jl)=zabscomp(jl) + & @@ -2302,6 +2306,7 @@ SUBROUTINE ham_rad_diag(kproma, kbdim, klev, krow, pxtm1) IF (nradang(1)/=0 .AND. nradang(2)/=0) THEN !>>SF #458 (replacing WHERE statements) + ! PLS - TODO: SAFE DIVISION ll1(1:kproma) = (tau_2d(nradang(1))%ptr(1:kproma,krow)>zeps) & .AND. (tau_2d(nradang(2))%ptr(1:kproma,krow)>zeps) diff --git a/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 b/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 index 6aaf5e1..43b518b 100755 --- a/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 +++ b/ifs-source/arpifs/m7/module/mo_tracer_processes.F90 @@ -139,7 +139,8 @@ SUBROUTINE xt_conv_massfix (kproma, kbdim, klev, & DO jk=1, klev DO jl=1, kproma - IF (ABS(zdxtdtsum(jl)) > zeps) THEN + ! PLS - TODO: SAFE DIVISION + IF (ABS(zdxtdtsum(jl)) > zeps) THEN zxttefix=-((ABS(zxtte(jl,jk)*zdpg(jl,jk))/zdxtdtsum(jl))*zdxtdt(jl))/zdpg(jl,jk) 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 fa227db..1cbdf4a 100644 --- a/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 +++ b/ifs-source/arpifs/m7/phys_ec/yoe_aer_activ.F90 @@ -520,6 +520,7 @@ SUBROUTINE AER_ACTIV_MORALES_NENES_FULL(KIDIA, KFDIA, KTDIA, KLON, KLEV, LCLOUD, NCL(JL) = NNACL(JL) NH2SO4(JL) = NSO4(JL) - NNA2SO4(JL) + ! PLS - TODO: SAFE DIVISION 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 f9e8b1207296dc95a583a2daea23ddba8588d739 Mon Sep 17 00:00:00 2001 From: Philippe Le Sager Date: Wed, 3 Jun 2026 15:03:56 +0200 Subject: [PATCH 5/6] Add M7 threshold to log --- ifs-source/arpifs/m7/hamm7_init.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ifs-source/arpifs/m7/hamm7_init.F90 b/ifs-source/arpifs/m7/hamm7_init.F90 index 0736abb..a6d0a67 100644 --- a/ifs-source/arpifs/m7/hamm7_init.F90 +++ b/ifs-source/arpifs/m7/hamm7_init.F90 @@ -55,6 +55,7 @@ SUBROUTINE hamm7_init(YGFL, YRRIP, CHEM_SCHEME) USE YOMLUN, ONLY : NULOUT ! --- M7 modules -------------------------------------------------------------- +USE MO_KIND, ONLY: THRESHOLD USE MO_TIME_CONTROL, ONLY: init_mo_time_control USE MO_SPECIES, ONLY: speclist ! tracer species in HAM USE MO_HAM, ONLY: & @@ -478,6 +479,8 @@ SUBROUTINE hamm7_init(YGFL, YRRIP, CHEM_SCHEME) ! -- LOG WRITE(NULOUT,'("====== HAMM7_INIT ===== ")') +WRITE(NULOUT,*) "TOLERANCE (EPS): ", THRESHOLD + WRITE(NULOUT,'("Number of size classes:", I3)') znclass WRITE(NULOUT,'(" class# / IFS id / HM7 id / IFSNAME / M7NAME ")') DO J_CLASS = 1,NCLASS From 1477b1b8fc83a1e14c3cb90d20b8c9600993aac2 Mon Sep 17 00:00:00 2001 From: Philippe Le Sager Date: Sun, 7 Jun 2026 23:44:26 +0200 Subject: [PATCH 6/6] Use erfc to compute CDF complement --- ifs-source/arpifs/m7/module/mo_ham_m7.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 index 7087165..574edcb 100644 --- a/ifs-source/arpifs/m7/module/mo_ham_m7.F90 +++ b/ifs-source/arpifs/m7/module/mo_ham_m7.F90 @@ -152,7 +152,7 @@ SUBROUTINE m7_cumulative_normal ( arg, presult, ccum ) REAL(dp), INTENT(OUT) :: presult presult = 0.5_dp * (1.0_dp + erf(arg / sqrt(2.0_dp))) - ccum = 1.0_dp - presult + ccum = 0.5_dp * erfc(arg / sqrt(2.0_dp)) END SUBROUTINE m7_cumulative_normal