Lima execution in mixed precision

UMR-CNRM · Nov 5, 2024 · 4cad834 · 4cad834
1 parent 9572d32
commit 4cad834
Show file tree

Hide file tree

Showing 8 changed files with 419 additions and 31 deletions.
diff --git a/build/with_fcm/arch/arch-gnu32.env b/build/with_fcm/arch/arch-gnu32.env
diff --git a/build/with_fcm/arch/arch-gnu32.fcm b/build/with_fcm/arch/arch-gnu32.fcm
@@ -0,0 +1,27 @@
+# Compilation
+$FCOMPILER     =     gfortran
+$BASE_FFLAGS   =     -c -fPIC -ffree-line-length-none -fopenmp -fconvert=swap -fallow-argument-mismatch
+$PROD_FFLAGS   =     -O2
+$DEV_FFLAGS    =     -O1
+$DEBUG_FFLAGS  =     -O0 -g -fbounds-check -finit-real=snan -ffpe-trap=invalid,zero,overflow
+$CCOMPILER     =     gcc
+$BASE_CFLAGS   =     -c -fPIC -fopenmp
+$PROD_CFLAGS   =     -O2
+$DEV_CFLAGS    =     -O1
+$DEBUG_CFLAGS  =     -fbounds-check
+$OMP_FFLAGS    =
+
+# Preprocessor
+$FPP_FLAGS     =     LINUX LITTLE_ENDIAN LITTLE REPRO48 PARKIND1_SINGLE
+$CPP_FLAGS     =     LINUX LITTLE_ENDIAN LITTLE PARKIND1_SINGLE
+$FPP_FLAGS_TESTPROGS = WITHOUT_CXXDEMANGLE USE_OPENMP
+
+# Linker
+$LINK          =     gfortran
+$BASE_LD       =     -fPIC  -fdefault-real-8 -fdefault-double-8 -fopenmp
+$OMP_LD        =
+$LD_EXE_TO_SHARED = -shared
+
+# Other
+$AR            =     ar
+
diff --git a/src/common/micro/lima.F90 b/src/common/micro/lima.F90
@@ -149,12 +149,12 @@ SUBROUTINE LIMA ( D, CST, ICED, ICEP, ELECD, ELECP, BUCONF, TBUDGETS, KBUDGETS,
 REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3))      :: ZCCT, ZCRT, ZCIT, ZCST, ZCGT, ZCHT
 REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3))      :: ZTHS, ZRVS, ZRCS, ZRRS, ZRIS, ZRSS, ZRGS, ZRHS
 REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3))      :: ZCCS, ZCRS, ZCIS, ZCSS, ZCGS, ZCHS
-REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NCCN) :: ZCCNFT, ZCCNAT
-REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NCCN) :: ZCCNFS, ZCCNAS
-REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NIFN) :: ZIFNFT, ZIFNNT
-REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NIFN) :: ZIFNFS, ZIFNNS
-REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NIMM) :: ZIMMNT
-REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NIMM) :: ZIMMNS
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NMOD_CCN) :: ZCCNFT, ZCCNAT
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NMOD_CCN) :: ZCCNFS, ZCCNAS
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NMOD_IFN) :: ZIFNFT, ZIFNNT
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NMOD_IFN) :: ZIFNFS, ZIFNNS
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NMOD_IMM) :: ZIMMNT
+REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3),NMOD_IMM) :: ZIMMNS
 REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3))      :: ZHOMFT
 REAL, DIMENSION(SIZE(PRT,1),SIZE(PRT,2),SIZE(PRT,3))      :: ZHOMFS
 

diff --git a/src/common/micro/mode_ini_lima_cold_mixed.F90 b/src/common/micro/mode_ini_lima_cold_mixed.F90
@@ -52,6 +52,7 @@ SUBROUTINE INI_LIMA_COLD_MIXED (PTSTEP, PDZMIN)
 USE MODD_PARAM_LIMA_MIXED
 !
 use mode_msg
+USE MODD_PRECISION, ONLY: MNHREAL64
 !
 USE MODE_LIMA_FUNCTIONS, ONLY: MOMG, GAUHER
 USE MODI_GAMMA
@@ -315,13 +316,13 @@ SUBROUTINE INI_LIMA_COLD_MIXED (PTSTEP, PDZMIN)
 !Cas GAMMAGEN
    XALPHAS = .214   ! Generalized gamma law
    XNUS    = 43.7   ! Generalized gamma law
-   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
-                            ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )
+   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(DBLE(XNUS + 2._MNHREAL64/XALPHAS))*GAMMA(DBLE(XNUS + 4._MNHREAL64/XALPHAS)) ) / &
+                            ( 8._MNHREAL64* GAMMA(DBLE(XNUS + 1._MNHREAL64/XALPHAS))*GAMMA(DBLE(XNUS + 3._MNHREAL64/XALPHAS)) ) )
 ELSE IF (NMOM_S.EQ.2) THEN
    XALPHAS = 1.0  ! Gamma law 
    XNUS    = 2.0  !
-   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
-                            ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )
+   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(DBLE(XNUS + 2._MNHREAL64/XALPHAS))*GAMMA(DBLE(XNUS + 4._MNHREAL64/XALPHAS)) ) / &
+                            ( 8._MNHREAL64* GAMMA(DBLE(XNUS + 1._MNHREAL64/XALPHAS))*GAMMA(DBLE(XNUS + 3._MNHREAL64/XALPHAS)) ) )
 ELSE
    XALPHAS = 1.0  ! Exponential law
    XNUS    = 1.0  ! Exponential law

diff --git a/src/common/micro/mode_ini_rain_ice.F90 b/src/common/micro/mode_ini_rain_ice.F90
@@ -99,6 +99,7 @@ SUBROUTINE INI_RAIN_ICE ( KLUOUT, PTSTEP, PDZMIN, KSPLITR, HCLOUD )
                         & XFRACM90, XFRMIN_NAM, XRDEPGRED_NAM, XRDEPSRED_NAM
 USE MODD_RAIN_ICE_DESCR_n
 USE MODD_RAIN_ICE_PARAM_n
+USE MODD_PRECISION, ONLY: MNHREAL64
 !
 USE MODI_GAMMA
 USE MODI_GAMMA_INC
@@ -392,8 +393,8 @@ SUBROUTINE INI_RAIN_ICE ( KLUOUT, PTSTEP, PDZMIN, KSPLITR, HCLOUD )
 !Cas GAMMAGEN
    XALPHAS = .214   ! Generalized gamma law
    XNUS    = 43.7   ! Generalized gamma law
-   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
-                            ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )
+   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(DBLE(XNUS + 2._MNHREAL64/XALPHAS))*GAMMA(DBLE(XNUS + 4._MNHREAL64/XALPHAS)) ) / &
+                            ( 8._MNHREAL64* GAMMA(DBLE(XNUS + 1._MNHREAL64/XALPHAS))*GAMMA(DBLE(XNUS + 3._MNHREAL64/XALPHAS)) ) )
 ELSE
    XALPHAS = 1.0  ! Exponential law
    XNUS    = 1.0  ! Exponential law

diff --git a/src/common/micro/mode_lima_init_ccn_activation_spectrum.F90 b/src/common/micro/mode_lima_init_ccn_activation_spectrum.F90
@@ -39,6 +39,8 @@ SUBROUTINE LIMA_INIT_CCN_ACTIVATION_SPECTRUM (CTYPE_CCN,XD,XSIGMA,XLIMIT_FACTOR,
 USE MODI_GAMMA_INC
 USE MODI_HYPGEO
 USE MODI_HYPSER
+USE MODD_PRECISION, ONLY: MNHREAL64
+USE MODI_MINPACK
 !
 IMPLICIT NONE
 !
@@ -55,13 +57,13 @@ SUBROUTINE LIMA_INIT_CCN_ACTIVATION_SPECTRUM (CTYPE_CCN,XD,XSIGMA,XLIMIT_FACTOR,
 !
 !*       0.2   Declarations of local variables :
 !
-INTEGER, PARAMETER            :: M = 1000        ! Number of points (S,Nccn) used to fit the spectra
-INTEGER, PARAMETER            :: N = 3          ! Number of parameters to adjust
-REAL, DIMENSION(N)            :: PARAMS         ! Parameters to adjust by the LM algorithm (k, mu, beta)
-REAL, DIMENSION(M)            :: FVEC           ! Array to store the distance between theoretical and fitted spectra
+INTEGER(KIND=4), PARAMETER    :: M = 1000        ! Number of points (S,Nccn) used to fit the spectra
+INTEGER(KIND=4), PARAMETER    :: N = 3          ! Number of parameters to adjust
+REAL(KIND=MNHREAL64), DIMENSION(N)            :: PARAMS         ! Parameters to adjust by the LM algorithm (k, mu, beta)
+REAL(KIND=MNHREAL64), DIMENSION(M)            :: FVEC           ! Array to store the distance between theoretical and fitted spectra
 INTEGER                       :: IFLAG          ! 
-INTEGER                       :: INFO           ! 
-REAL                          :: TOL = 1.E-16   ! Fit precision required
+INTEGER(KIND=4)               :: INFO           ! 
+REAL(KIND=MNHREAL64)          :: TOL = 1.E-16   ! Fit precision required
 !
 INTEGER                       :: II, IJ         ! Loop indices
 !
@@ -152,8 +154,9 @@ SUBROUTINE LIMA_INIT_CCN_ACTIVATION_SPECTRUM (CTYPE_CCN,XD,XSIGMA,XLIMIT_FACTOR,
 !*       3.     Compute C, k, mu, beta, using the Levenberg-Marquardt algorithm
 !	        ---------------------------------------------------------------
 !
-   PARAMS(1:3) = (/ 1., 1., 1000. /)
+   PARAMS(1:3) = (/ 1._MNHREAL64, 1._MNHREAL64, 1000._MNHREAL64 /)
    IFLAG = 1
+   !lmdif1 uses KIND 8 reals and KIND 4 integers
    call lmdif1 ( DISTANCE, M, N, PARAMS, FVEC, TOL, INFO )
 !
    XLIMIT_FACTOR = gamma(PARAMS(2))*PARAMS(3)**(PARAMS(1)/2)/gamma(1+PARAMS(1)/2)/gamma(PARAMS(2)-PARAMS(1)/2)
@@ -333,6 +336,7 @@ FUNCTION DSDD(XD,XDDRY,XKAPPA, XT)  RESULT(DS)
 !*       0. DECLARATIONS
 !
     USE MODD_CST, ONLY : XMV, XAVOGADRO, XBOLTZ, XRHOLW
+    USE MODD_PRECISION, ONLY: MNHREAL64
 !
     IMPLICIT NONE
 !
@@ -348,10 +352,12 @@ FUNCTION DSDD(XD,XDDRY,XKAPPA, XT)  RESULT(DS)
 !*       0.2 declarations of local variables
 !
     REAL              :: XA     ! factor inside the exponential
+    REAL(KIND=MNHREAL64)      :: Z
 !    
     XA = 4 * 0.072 * XMV / XAVOGADRO / XBOLTZ / XT / XRHOLW
-    DS = (XD**3-XDDRY**3) * (XD**3-(1-XKAPPA)*XDDRY**3) * XA - 3. * XKAPPA * XD**4 * XDDRY**3
-    DS = DS * EXP(XA/XD) / (XD**3-(1-XKAPPA)*XDDRY**3)**2
+    Z = (XD**3-XDDRY**3) * (XD**3-(1._MNHREAL64-XKAPPA)*XDDRY**3) * XA - 3._MNHREAL64 * XKAPPA * DBLE(XD)**4 * DBLE(XDDRY)**3
+    Z = Z * EXP(XA/XD) / (XD**3-(1-XKAPPA)*XDDRY**3)**2
+    DS = Z
 !
 END FUNCTION DSDD
 !
@@ -396,11 +402,12 @@ SUBROUTINE DISTANCE(M,N,X,FVEC,IFLAG)
 !
 !*       0.1 declarations of arguments and result
 !
-    integer, intent(in) :: M
-    integer, intent(in) :: N
-    real, intent(in) ::    X(N)
-    real, intent(out) ::    FVEC(M)
-    integer, intent(inout) :: IFLAG
+    !DISTANCE must use KIND 8 reals and KIND 4 integers to be used by LMDIF1
+    integer(KIND=4), intent(in) :: M
+    integer(KIND=4), intent(in) :: N
+    real(KIND=MNHREAL64), intent(in) ::    X(N)
+    real(KIND=MNHREAL64), intent(out) ::    FVEC(M)
+    integer(KIND=4), intent(inout) :: IFLAG
 !
 !*       0.2 declarations of local variables
 !
@@ -416,7 +423,7 @@ SUBROUTINE DISTANCE(M,N,X,FVEC,IFLAG)
        DO I=1, M
           ! XS in "no units", ie XS=0.01 for a 1% suersaturation
           !       ZW= C * (XS(I)/100)**X(1) * HYPGEO(X(2),X(1)/2,X(1)/2+1,X(3),XS(I)/100)
-          ZW= C * (XS(I))**X(1) * HYPGEO(X(2),X(1)/2,X(1)/2+1,X(3),XS(I))
+          ZW= C * (XS(I))**X(1) * HYPGEO(REAL(X(2)), REAL(X(1)/2), REAL(X(1)/2+1), REAL(X(3)), REAL(XS(I)))
 !!$       IF (X(3)*(XS(I)/100)**2 .LT. 0.98) THEN
 !!$          CALL HYPSER(X(2),X(1)/2,X(1)/2+1,-X(3)*(XS(I)/100)**2,ZW2)
 !!$          print *, "args= ", X(2), X(1)/2, X(1)/2+1, -X(3)*(XS(I)/100)**2, " hypser = ", ZW2