From c2f58d86443006eaefa8ed8176dbd9f2f4347001 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Vi=C3=A9?= <benoit.vie@meteo.fr>
Date: Tue, 24 May 2022 19:33:32 +0200
Subject: [PATCH] bugfix

---
 src/MNH/ice4_fast_rs.f90                 |   4 +-
 src/MNH/ice4_slow.f90                    |   2 +-
 src/MNH/lima_cold_slow_processes.f90     |   2 +-
 src/MNH/lima_conversion_melting_snow.f90 |   4 +-
 src/MNH/lima_ice_snow_deposition.f90     | 233 -----------------------
 src/MNH/lima_mixed_fast_processes.f90    |   4 +-
 src/MNH/lima_sedimentation.f90           |   2 +-
 src/MNH/lima_snow_deposition.f90         |   6 +-
 8 files changed, 12 insertions(+), 245 deletions(-)
 delete mode 100644 src/MNH/lima_ice_snow_deposition.f90

diff --git a/src/MNH/ice4_fast_rs.f90 b/src/MNH/ice4_fast_rs.f90
index 11dbd67f5..d44b1bcc7 100644
--- a/src/MNH/ice4_fast_rs.f90
+++ b/src/MNH/ice4_fast_rs.f90
@@ -173,7 +173,7 @@ ELSE
                            *(XESTT-PRS_TEND(:, IFREEZ1))/(XRV*PT(:))           )
     PRS_TEND(:, IFREEZ1)=PRS_TEND(:, IFREEZ1)* PRST(:) * ( X0DEPS*       PLBDAS(:)**XEX0DEPS +     &
          X1DEPS*PCJ(:)*PLBDAS(:) **(XBS+XEX1DEPS)* &
-         (1+(XFVELOS/(2*PLBDAS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS))/ &
+         (1+0.5*(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS))/ &
          ( PRHODREF(:)*(XLMTT-XCL*(XTT-PT(:))) )
     PRS_TEND(:, IFREEZ2)=(PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PT(:)))   ) / &
                           ( PRHODREF(:)*(XLMTT-XCL*(XTT-PT(:))) )
@@ -503,7 +503,7 @@ ELSE
     PRSMLTG(:)  = XFSCVMG*MAX( 0.0,( -PRSMLTG(:) *             &
          PRST(:)*PRHODREF(:) *    &
          ( X0DEPS       *PLBDAS(:)**XEX0DEPS +     &
-         X1DEPS*PCJ(:)*(1+(XFVELOS/(2*PLBDAS(:))**XALPHAS))**(XNUS+XEX1DEPS/XALPHAS)*((PLBDAS(:))**(XBS+XEX1DEPS))) -   &
+         X1DEPS*PCJ(:)*(1+0.5*(XFVELOS/PLBDAS(:))**XALPHAS)**(XNUS+XEX1DEPS/XALPHAS)*PLBDAS(:)**(XBS+XEX1DEPS)) -   &
          ( PRS_TEND(:, IRCRIMS) + PRS_TEND(:, IRRACCS)) *       &
          ( PRHODREF(:)*XCL*(XTT-PT(:))) ) /    &
          ( PRHODREF(:)*XLMTT ) )
diff --git a/src/MNH/ice4_slow.f90 b/src/MNH/ice4_slow.f90
index b9e202b14..b9325050f 100644
--- a/src/MNH/ice4_slow.f90
+++ b/src/MNH/ice4_slow.f90
@@ -176,7 +176,7 @@ ELSE
   DO JL=1, KSIZE
      IF (ZMASK(JL)==1.) THEN
         PRVDEPS(JL) = ( PRST(JL)*PSSI(JL)/PAI(JL) ) * &
-             ( X0DEPS*PLBDAS(JL)**XEX0DEPS + X1DEPS*PCJ(JL) * (1+(XFVELOS/(2*PLBDAS(JL)))**XALPHAS)**(-XNUS+XEX1DEPS) &
+             ( X0DEPS*PLBDAS(JL)**XEX0DEPS + X1DEPS*PCJ(JL) * (1+0.5*(XFVELOS/PLBDAS(JL))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) &
              *(PLBDAS(JL))**(XBS+XEX1DEPS) )
      END IF
   END DO
diff --git a/src/MNH/lima_cold_slow_processes.f90 b/src/MNH/lima_cold_slow_processes.f90
index a74b94686..29c1aecbe 100644
--- a/src/MNH/lima_cold_slow_processes.f90
+++ b/src/MNH/lima_cold_slow_processes.f90
@@ -393,7 +393,7 @@ IF( IMICRO >= 1 ) THEN
       WHERE ( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>ZRTMIN(5)) )
          ZZW(:) = ( ZRST(:)*ZSSI(:)/(ZAI(:)) ) *                                            &
               ( X0DEPS*ZLBDAS(:)**XEX0DEPS +                                                &
-              (X1DEPS*ZCJ(:)*(1+(XFVELOS/(2.*ZLBDAS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) * &
+              (X1DEPS*ZCJ(:)*(1+0.5*(XFVELOS/ZLBDAS(:))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) * &
               (ZLBDAS(:))**(XEX1DEPS+XBS)))
 
          ZZW(:) =    MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
diff --git a/src/MNH/lima_conversion_melting_snow.f90 b/src/MNH/lima_conversion_melting_snow.f90
index 307db0255..d1aa60aa4 100644
--- a/src/MNH/lima_conversion_melting_snow.f90
+++ b/src/MNH/lima_conversion_melting_snow.f90
@@ -107,10 +107,10 @@ WHERE( (PRST(:)>XRTMIN(5)) .AND. (PT(:)>XTT) .AND. LDCOMPUTE(:) )
 !
 ! compute RSMLT
 !
-   ZW(:)  = XFSCVMG*MAX( 0.0,( -ZW(:) * PRHODREF(:) * PRST(:) *          &
+   ZW(:)  = XFSCVMG*MAX( 0.0,( -ZW(:) * PRST(:) *                        &
                                ( X0DEPS*PLBDS(:)**XEX0DEPS +             &
                                  X1DEPS*PCJ(:)*PLBDS(:)**(XEX1DEPS+XBS)* &
-                                   (1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) ))
+                                   (1+0.5*(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) ))
 ! On ne tient pas compte de la collection de pluie et gouttelettes par la neige si T>0 !!!! 
 ! Note that no heat is exchanged because the graupeln produced are still icy!!!
    P_RS_CMEL(:) = - ZW(:)
diff --git a/src/MNH/lima_ice_snow_deposition.f90 b/src/MNH/lima_ice_snow_deposition.f90
deleted file mode 100644
index 4968a192d..000000000
--- a/src/MNH/lima_ice_snow_deposition.f90
+++ /dev/null
@@ -1,233 +0,0 @@
-!MNH_LIC Copyright 2013-2018 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
-!MNH_LIC for details. version 1.
-!      #####################
-       MODULE MODI_LIMA_ICE_SNOW_DEPOSITION
-!      #####################
-!
-INTERFACE
-      SUBROUTINE LIMA_ICE_SNOW_DEPOSITION (PTSTEP, LDCOMPUTE,                 &
-                                           PRHODREF, PSSI, PAI, PCJ, PLSFACT, &
-                                           PRIT, PRST, PCIT, PLBDI, PLBDS,    &
-                                           P_RI_CNVI, P_CI_CNVI,              &
-                                           P_TH_DEPS, P_RS_DEPS,              &
-                                           P_RI_CNVS, P_CI_CNVS,              &
-                                           PA_TH, PA_RV, PA_RI, PA_CI, PA_RS  )
-!
-REAL,                 INTENT(IN)    :: PTSTEP
-LOGICAL, DIMENSION(:),INTENT(IN)    :: LDCOMPUTE
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PRHODREF! Reference density
-REAL, DIMENSION(:),   INTENT(IN)    :: PSSI  ! abs. pressure at time t
-REAL, DIMENSION(:),   INTENT(IN)    :: PAI  ! abs. pressure at time t
-REAL, DIMENSION(:),   INTENT(IN)    :: PCJ  ! abs. pressure at time t
-REAL, DIMENSION(:),   INTENT(IN)    :: PLSFACT  ! abs. pressure at time t
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PRIT    ! Cloud ice m.r. at t 
-REAL, DIMENSION(:),   INTENT(IN)    :: PRST    ! Snow/aggregate m.r. at t 
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PCIT    ! Ice crystal C. at t
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PLBDI    ! Graupel m.r. at t 
-REAL, DIMENSION(:),   INTENT(IN)    :: PLBDS    ! Graupel m.r. at t 
-!
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_RI_CNVI
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_CI_CNVI
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_TH_DEPS
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_RS_DEPS
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_RI_CNVS
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_CI_CNVS
-!
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_TH
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_RV
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_RI
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_CI
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_RS
-!
-END SUBROUTINE LIMA_ICE_SNOW_DEPOSITION
-END INTERFACE
-END MODULE MODI_LIMA_ICE_SNOW_DEPOSITION
-!
-!     ##########################################################################
-SUBROUTINE LIMA_ICE_SNOW_DEPOSITION (PTSTEP, LDCOMPUTE,                        &
-                                     PRHODREF,  PSSI, PAI, PCJ, PLSFACT,       &
-                                     PRIT, PRST, PCIT, PLBDI, PLBDS,           &
-                                     P_RI_CNVI, P_CI_CNVI,                     &
-                                     P_TH_DEPS, P_RS_DEPS,                     &
-                                     P_RI_CNVS, P_CI_CNVS,                     &
-                                     PA_TH, PA_RV, PA_RI, PA_CI, PA_RS         )
-!     ##########################################################################
-!
-!!    PURPOSE
-!!    -------
-!!      The purpose of this routine is to compute the microphysical sources
-!!    for slow cold processes :
-!!      - conversion of snow to ice
-!!      - deposition of vapor on snow
-!!      - conversion of ice to snow (Harrington 1995)
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!      J.-M. Cohard     * Laboratoire d'Aerologie*
-!!      J.-P. Pinty      * Laboratoire d'Aerologie*
-!!      S.    Berthet    * Laboratoire d'Aerologie*
-!!      B.    Vié        * CNRM *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original             15/03/2018
-!  J. Wurtz       03/2022: new snow characteristics
-!!
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-USE MODD_PARAM_LIMA,      ONLY : XRTMIN, XCTMIN, XALPHAI, XALPHAS, XNUI, XNUS
-USE MODD_PARAM_LIMA_COLD, ONLY : XCXS, XNS, XBS, &
-                                 XLBDAS_MAX, XDSCNVI_LIM, XLBDASCNVI_MAX,     &
-                                 XC0DEPSI, XC1DEPSI, XR0DEPSI, XR1DEPSI,      &
-                                 XSCFAC, X1DEPS, X0DEPS, XEX1DEPS, XEX0DEPS,  &
-                                 XDICNVS_LIM, XLBDAICNVS_LIM,                 &
-                                 XC0DEPIS, XC1DEPIS, XR0DEPIS, XR1DEPIS,      &
-                                 XCOLEXIS, XAGGS_CLARGE1, XAGGS_CLARGE2,      &
-                                 XAGGS_RLARGE1, XAGGS_RLARGE2, XFVELOS
-
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of dummy arguments :
-!
-REAL,                 INTENT(IN)    :: PTSTEP
-LOGICAL, DIMENSION(:),INTENT(IN)    :: LDCOMPUTE
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PRHODREF! Reference density
-REAL, DIMENSION(:),   INTENT(IN)    :: PSSI  ! abs. pressure at time t
-REAL, DIMENSION(:),   INTENT(IN)    :: PAI  ! abs. pressure at time t
-REAL, DIMENSION(:),   INTENT(IN)    :: PCJ  ! abs. pressure at time t
-REAL, DIMENSION(:),   INTENT(IN)    :: PLSFACT  ! abs. pressure at time t
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PRIT    ! Cloud ice m.r. at t 
-REAL, DIMENSION(:),   INTENT(IN)    :: PRST    ! Snow/aggregate m.r. at t 
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PCIT    ! Ice crystal C. at t
-!
-REAL, DIMENSION(:),   INTENT(IN)    :: PLBDI    ! Graupel m.r. at t 
-REAL, DIMENSION(:),   INTENT(IN)    :: PLBDS    ! Graupel m.r. at t 
-!
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_RI_CNVI
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_CI_CNVI
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_TH_DEPS
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_RS_DEPS
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_RI_CNVS
-REAL, DIMENSION(:),   INTENT(INOUT) :: P_CI_CNVS
-!
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_TH
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_RV
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_RI
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_CI
-REAL, DIMENSION(:),   INTENT(INOUT) :: PA_RS
-!
-!*       0.2   Declarations of local variables :
-!
-LOGICAL, DIMENSION(SIZE(PRHODREF)) :: GMICRO ! Computations only where necessary
-REAL,    DIMENSION(SIZE(PRHODREF)) :: ZZW, ZZW2, ZZX ! Work array
-!
-!
-!-------------------------------------------------------------------------------
-!
-P_RI_CNVI(:) = 0.
-P_CI_CNVI(:) = 0.
-P_TH_DEPS(:) = 0.
-P_RS_DEPS(:) = 0.
-P_RI_CNVS(:) = 0.
-P_CI_CNVS(:) = 0.
-!
-! Physical limitations
-!
-!
-! Looking for regions where computations are necessary
-!
-GMICRO(:) = .FALSE.
-GMICRO(:) = LDCOMPUTE(:) .AND. &
-     (PRIT(:)>XRTMIN(4)      .OR.  &
-      PRST(:)>XRTMIN(5))
-!
-!
-WHERE( GMICRO )
-!
-!*       2.1    Conversion of snow to r_i: RSCNVI
-!        ----------------------------------------
-!
-!
-   ZZW2(:) = 0.0
-   ZZW(:) = 0.0
-   WHERE ( PLBDS(:)<XLBDASCNVI_MAX .AND. (PRST(:)>XRTMIN(5)) &
-                                   .AND. (PSSI(:)<0.0)       )
-      ZZW(:) = (PLBDS(:)*XDSCNVI_LIM)**(XALPHAS)
-      ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XNS*PRST(:)*PLBDS(:)**XBS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
-!
-      ZZW(:) = ( XR0DEPSI+XR1DEPSI*PCJ(:) )*ZZX(:)
-!
-      ZZW2(:) = ZZW(:)*( XC0DEPSI+XC1DEPSI*PCJ(:) )/( XR0DEPSI+XR1DEPSI*PCJ(:) )
-   END WHERE
-!
-   P_RI_CNVI(:) = ZZW(:)
-   P_CI_CNVI(:) = ZZW2(:)
-!
-   PA_RI(:) = PA_RI(:) + P_RI_CNVI(:)
-   PA_CI(:) = PA_CI(:) + P_CI_CNVI(:)
-   PA_RS(:) = PA_RS(:) - P_RI_CNVI(:)
-!
-!
-!*       2.2    Deposition of water vapor on r_s: RVDEPS
-!        -----------------------------------------------
-!
-!
-   ZZW(:) = 0.0
-   WHERE ( (PRST(:)>XRTMIN(5)) )
-      ZZW(:) =( PRST(:)*PSSI(:)/PAI(:) ) *                               &
-           ( X0DEPS*PLBDS(:)**XEX0DEPS +                &
-           ( X1DEPS*PCJ(:)*(PLBDS(:))**(XBS+XEX1DEPS) * &
-                    (1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)))
-      ZZW(:) =    ZZW(:)*(0.5+SIGN(0.5,ZZW(:))) - ABS(ZZW(:))*(0.5-SIGN(0.5,ZZW(:)))
-   END WHERE
-!
-   P_RS_DEPS(:) = ZZW(:)
-   P_TH_DEPS(:) = P_RS_DEPS(:) * PLSFACT(:)
-!
-   PA_TH(:) = PA_TH(:) + P_TH_DEPS(:)
-   PA_RV(:) = PA_RV(:) - P_RS_DEPS(:) 
-   PA_RS(:) = PA_RS(:) + P_RS_DEPS(:) 
-!
-!
-!*       2.3    Conversion of pristine ice to r_s: RICNVS
-!        ------------------------------------------------
-!
-!
-   ZZW(:) = 0.0
-   ZZW2(:) = 0.0
-   WHERE ( (PLBDI(:)<XLBDAICNVS_LIM) .AND. (PCIT(:)>XCTMIN(4)) &
-                                     .AND. (PSSI(:)>0.0)       )
-      ZZW(:) = (PLBDI(:)*XDICNVS_LIM)**(XALPHAI)
-      ZZX(:) = ( PSSI(:)/PAI(:) )*PCIT(:) * (ZZW(:)**XNUI) *EXP(-ZZW(:))
-!
-      ZZW(:) = ( XR0DEPIS + XR1DEPIS*PCJ(:) )*ZZX(:)                             
-!
-      ZZW2(:) = ZZW(:) * (XC0DEPIS+XC1DEPIS*PCJ(:)) / (XR0DEPIS+XR1DEPIS*PCJ(:))
-   END WHERE
-!
-P_RI_CNVS(:) = - ZZW(:)
-P_CI_CNVS(:) = - ZZW2(:)
-!
-PA_RI(:) = PA_RI(:) + P_RI_CNVS(:)
-PA_CI(:) = PA_CI(:) + P_CI_CNVS(:)
-PA_RS(:) = PA_RS(:) - P_RI_CNVS(:)
-!
-!
-END WHERE
-!
-!
-END SUBROUTINE LIMA_ICE_SNOW_DEPOSITION
diff --git a/src/MNH/lima_mixed_fast_processes.f90 b/src/MNH/lima_mixed_fast_processes.f90
index 94b622ef3..056386015 100644
--- a/src/MNH/lima_mixed_fast_processes.f90
+++ b/src/MNH/lima_mixed_fast_processes.f90
@@ -923,8 +923,8 @@ WHERE( (PRST1D(:)>XRTMIN(5)) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP) .AND. (PZT(:)>XT
 ! compute RSMLT
 !
   ZZW(:)  = MIN( PRSS1D(:), XFSCVMG*MAX( 0.0,( -ZZW(:) *             &
-                          PRHODREF(:) * PRST1D(:)*( X0DEPS*            PLBDAS(:)**XEX0DEPS +     &
-                          X1DEPS*PCJ(:)*(1+(XFVELOS/(2.*PLBDAS(:)))**XALPHAS) &
+                          PRST1D(:)*( X0DEPS*            PLBDAS(:)**XEX0DEPS +     &
+                          X1DEPS*PCJ(:)*(1+0.5*(XFVELOS/PLBDAS(:))**XALPHAS) &
                           **(-XNUS+XEX1DEPS/XALPHAS)*(PLBDAS(:))**(XEX1DEPS+XBS))-    &
                                     ( ZZW1(:,1)+ZZW1(:,4) ) *       &
                              ( PRHODREF(:)*XCL*(XTT-PZT(:))) ) /    &
diff --git a/src/MNH/lima_sedimentation.f90 b/src/MNH/lima_sedimentation.f90
index 8d48b776d..cd90504a1 100644
--- a/src/MNH/lima_sedimentation.f90
+++ b/src/MNH/lima_sedimentation.f90
@@ -195,7 +195,7 @@ DO JN = 1 ,  NSPLITSED(KID)
          END WHERE
          ZLBDA(:) = ZLBDA(:)*XTRANS_MP_GAMMAS
          ZZW(:) = XFSEDR(KID) * ZRHODREF(:)**(1.-XCEXVT)*ZRS(:)* &
-              (1 + (XFVELOS/ZLBDA(:))**XALPHAS)**(-XNUS+XEXSEDS/XALPHAS) * ZLBDA(:)**(XBS+XEXSEDS)
+              (1 + (XFVELOS/ZLBDA(:))**XALPHAS)**(-XNUS-(XD(KID)+XBS)/XALPHAS) * ZLBDA(:)**(-XD(KID))
       ELSE
          IF (KMOMENTS==1) ZLBDA(:) = XLB(KID) * ( ZRHODREF(:) * ZRS(:) )**XLBEX(KID)
          IF (KMOMENTS==2) ZLBDA(:) = ( XLB(KID)*ZCS(:) / ZRS(:) )**XLBEX(KID)
diff --git a/src/MNH/lima_snow_deposition.f90 b/src/MNH/lima_snow_deposition.f90
index fd28d8c6d..5a6ec35ca 100644
--- a/src/MNH/lima_snow_deposition.f90
+++ b/src/MNH/lima_snow_deposition.f90
@@ -125,10 +125,10 @@ IF (NMOM_I.EQ.1) THEN
 !
       ZZW(:) = 0.0
       WHERE ( (PRST(:)>XRTMIN(5)) )
-         ZZW(:) = PRHODREF(:) * PRST(:) * PSSI(:) / PAI(:) * &
+         ZZW(:) = PRST(:) * PSSI(:) / PAI(:) * &
               ( X0DEPS*PLBDS(:)**XEX0DEPS +                  &
                 X1DEPS*PLBDS(:)**(XEX1DEPS+XBS)*PCJ(:) *     &
-                     (1+(XFVELOS/(2.*PLBDS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) )
+                     (1+0.5*(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) )
          ZZW(:) =    ZZW(:)*(0.5+SIGN(0.5,ZZW(:))) - ABS(ZZW(:))*(0.5-SIGN(0.5,ZZW(:)))
       END WHERE
       P_RS_DEPS(:) = ZZW(:)
@@ -165,7 +165,7 @@ ELSE
          ZZW(:) = ( PRST(:)*PSSI(:)/(PAI(:)) ) *           &
               ( X0DEPS*PLBDS(:)**XEX0DEPS +                &
               ( X1DEPS*PCJ(:)*(PLBDS(:))**(XBS+XEX1DEPS) * &
-                   (1+(XFVELOS/(2.*PLBDS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) )
+                   (1+0.5*(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) )
          ZZW(:) =    ZZW(:)*(0.5+SIGN(0.5,ZZW(:))) - ABS(ZZW(:))*(0.5-SIGN(0.5,ZZW(:)))
       END WHERE
 !
-- 
GitLab