diff --git a/src/MNH/rain_ice_fast_rg.f90 b/src/MNH/rain_ice_fast_rg.f90
index 41ce09cbf7280e98c8b6095f2feddc1137f51027..7c2b0aeb7a600ad27138e4e99475ae6029b53dd8 100644
--- a/src/MNH/rain_ice_fast_rg.f90
+++ b/src/MNH/rain_ice_fast_rg.f90
@@ -7,6 +7,7 @@
! P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
! P. Wautelet 03/06/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+! P. Wautelet 05/06/2019: optimisations
!-----------------------------------------------------------------
MODULE MODE_RAIN_ICE_FAST_RG
@@ -78,19 +79,18 @@ REAL, DIMENSION(:), intent(out) :: PRWETG ! Wet growth rate of the g
!* 0.2 declaration of local variables
!
INTEGER :: IGDRY
-INTEGER :: JJ
+INTEGER :: JJ, JL
INTEGER, DIMENSION(size(PRHODREF)) :: I1
INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC1, IVEC2 ! Vectors of indices for interpolations
-LOGICAL, DIMENSION(size(PRHODREF)) :: GDRY ! Test where to compute dry growth
REAL, DIMENSION(size(PRHODREF)) :: ZZW ! Work array
REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1,ZVEC2,ZVEC3 ! Work vectors for interpolations
+REAL, DIMENSION(:), ALLOCATABLE :: ZVECLBDAG, ZVECLBDAR, ZVECLBDAS
REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!-------------------------------------------------------------------------------
!
!* 6.1 rain contact freezing
!
- ZZW1(:,3:4) = 0.0
WHERE( (PRIT(:)>XRTMIN(4)) .AND. (PRRT(:)>XRTMIN(3)) .AND. &
(PRIS(:)>0.0) .AND. (PRRS(:)>0.0) )
ZZW1(:,3) = MIN( PRIS(:),XICFRR * PRIT(:) & ! RICFRRG
@@ -120,11 +120,11 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!* 6.2 compute the Dry growth case
!
ZZW1(:,:) = 0.0
- WHERE( (PRGT(:)>XRTMIN(6)) .AND. ((PRCT(:)>XRTMIN(2) .AND. PRCS(:)>0.0)) )
+ WHERE( PRGT(:)>XRTMIN(6) .AND. PRCT(:)>XRTMIN(2) .AND. PRCS(:)>0.0 )
ZZW(:) = PLBDAG(:)**(XCXG-XDG-2.0) * PRHODREF(:)**(-XCEXVT)
ZZW1(:,1) = MIN( PRCS(:),XFCDRYG * PRCT(:) * ZZW(:) ) ! RCDRYG
END WHERE
- WHERE( (PRGT(:)>XRTMIN(6)) .AND. ((PRIT(:)>XRTMIN(4) .AND. PRIS(:)>0.0)) )
+ WHERE( (PRGT(:)>XRTMIN(6)) .AND. PRIT(:)>XRTMIN(4) .AND. PRIS(:)>0.0 )
ZZW(:) = PLBDAG(:)**(XCXG-XDG-2.0) * PRHODREF(:)**(-XCEXVT)
ZZW1(:,2) = MIN( PRIS(:),XFIDRYG * EXP( XCOLEXIG*(PZT(:)-XTT) ) &
* PRIT(:) * ZZW(:) ) ! RIDRYG
@@ -133,13 +133,10 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!* 6.2.1 accretion of aggregates on the graupeln
!
IGDRY = 0
- DO JJ = 1, SIZE(GDRY)
+ DO JJ = 1, SIZE(PRST)
IF ( PRST(JJ)>XRTMIN(5) .AND. PRGT(JJ)>XRTMIN(6) .AND. PRSS(JJ)>0.0 ) THEN
IGDRY = IGDRY + 1
I1(IGDRY) = JJ
- GDRY(JJ) = .TRUE.
- ELSE
- GDRY(JJ) = .FALSE.
END IF
END DO
@@ -147,6 +144,8 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!* 6.2.2 allocations
!
+ ALLOCATE(ZVECLBDAG(IGDRY))
+ ALLOCATE(ZVECLBDAS(IGDRY))
ALLOCATE(ZVEC1(IGDRY))
ALLOCATE(ZVEC2(IGDRY))
ALLOCATE(ZVEC3(IGDRY))
@@ -155,22 +154,20 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!* 6.2.3 select the (PLBDAG,PLBDAS) couplet
!
- DO JJ = 1, IGDRY
- ZVEC1(JJ) = PLBDAG(I1(JJ))
- ZVEC2(JJ) = PLBDAS(I1(JJ))
- END DO
+ ZVECLBDAG(1:IGDRY) = PLBDAG(I1(1:IGDRY))
+ ZVECLBDAS(1:IGDRY) = PLBDAS(I1(1:IGDRY))
!
!* 6.2.4 find the next lower indice for the PLBDAG and for the PLBDAS
! in the geometrical set of (Lbda_g,Lbda_s) couplet use to
! tabulate the SDRYG-kernel
!
ZVEC1(1:IGDRY) = MAX( 1.00001, MIN( REAL(NDRYLBDAG)-0.00001, &
- XDRYINTP1G * LOG( ZVEC1(1:IGDRY) ) + XDRYINTP2G ) )
+ XDRYINTP1G * LOG( ZVECLBDAG(1:IGDRY) ) + XDRYINTP2G ) )
IVEC1(1:IGDRY) = INT( ZVEC1(1:IGDRY) )
ZVEC1(1:IGDRY) = ZVEC1(1:IGDRY) - REAL( IVEC1(1:IGDRY) )
!
ZVEC2(1:IGDRY) = MAX( 1.00001, MIN( REAL(NDRYLBDAS)-0.00001, &
- XDRYINTP1S * LOG( ZVEC2(1:IGDRY) ) + XDRYINTP2S ) )
+ XDRYINTP1S * LOG( ZVECLBDAS(1:IGDRY) ) + XDRYINTP2S ) )
IVEC2(1:IGDRY) = INT( ZVEC2(1:IGDRY) )
ZVEC2(1:IGDRY) = ZVEC2(1:IGDRY) - REAL( IVEC2(1:IGDRY) )
!
@@ -185,20 +182,19 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
- XKER_SDRYG(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) &
* (ZVEC1(JJ) - 1.0)
END DO
- ZZW(:) = 0.
+!
DO JJ = 1, IGDRY
- ZZW(I1(JJ)) = ZVEC3(JJ)
+ JL = I1(JJ)
+ ZZW1(JL,3) = MIN( PRSS(JL),XFSDRYG*ZVEC3(JJ) & ! RSDRYG
+ * EXP( XCOLEXSG*(PZT(JL)-XTT) ) &
+ *( ZVECLBDAS(JJ)**(XCXS-XBS) )*( ZVECLBDAG(JJ)**XCXG ) &
+ *( PRHODREF(JL)**(-XCEXVT-1.) ) &
+ *( XLBSDRYG1/( ZVECLBDAG(JJ)**2 ) + &
+ XLBSDRYG2/( ZVECLBDAG(JJ) * ZVECLBDAS(JJ) ) + &
+ XLBSDRYG3/( ZVECLBDAS(JJ)**2) ) )
END DO
-!
- WHERE( GDRY(:) )
- ZZW1(:,3) = MIN( PRSS(:),XFSDRYG*ZZW(:) & ! RSDRYG
- * EXP( XCOLEXSG*(PZT(:)-XTT) ) &
- *( PLBDAS(:)**(XCXS-XBS) )*( PLBDAG(:)**XCXG ) &
- *( PRHODREF(:)**(-XCEXVT-1.) ) &
- *( XLBSDRYG1/( PLBDAG(:)**2 ) + &
- XLBSDRYG2/( PLBDAG(:) * PLBDAS(:) ) + &
- XLBSDRYG3/( PLBDAS(:)**2) ) )
- END WHERE
+ DEALLOCATE(ZVECLBDAS)
+ DEALLOCATE(ZVECLBDAG)
DEALLOCATE(IVEC2)
DEALLOCATE(IVEC1)
DEALLOCATE(ZVEC3)
@@ -209,13 +205,10 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!* 6.2.6 accretion of raindrops on the graupeln
!
IGDRY = 0
- DO JJ = 1, SIZE(GDRY)
+ DO JJ = 1, SIZE(PRRT)
IF ( PRRT(JJ)>XRTMIN(3) .AND. PRGT(JJ)>XRTMIN(6) .AND. PRRS(JJ)>0.0 ) THEN
IGDRY = IGDRY + 1
I1(IGDRY) = JJ
- GDRY(JJ) = .TRUE.
- ELSE
- GDRY(JJ) = .FALSE.
END IF
END DO
!
@@ -223,6 +216,8 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!* 6.2.7 allocations
!
+ ALLOCATE(ZVECLBDAG(IGDRY))
+ ALLOCATE(ZVECLBDAR(IGDRY))
ALLOCATE(ZVEC1(IGDRY))
ALLOCATE(ZVEC2(IGDRY))
ALLOCATE(ZVEC3(IGDRY))
@@ -231,22 +226,20 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!* 6.2.8 select the (PLBDAG,PLBDAR) couplet
!
- DO JJ = 1, IGDRY
- ZVEC1(JJ) = PLBDAG(I1(JJ))
- ZVEC2(JJ) = PLBDAR(I1(JJ))
- END DO
+ ZVECLBDAG(1:IGDRY) = PLBDAG(I1(1:IGDRY))
+ ZVECLBDAR(1:IGDRY) = PLBDAR(I1(1:IGDRY))
!
!* 6.2.9 find the next lower indice for the PLBDAG and for the PLBDAR
! in the geometrical set of (Lbda_g,Lbda_r) couplet use to
! tabulate the RDRYG-kernel
!
ZVEC1(1:IGDRY) = MAX( 1.00001, MIN( REAL(NDRYLBDAG)-0.00001, &
- XDRYINTP1G * LOG( ZVEC1(1:IGDRY) ) + XDRYINTP2G ) )
+ XDRYINTP1G * LOG( ZVECLBDAG(1:IGDRY) ) + XDRYINTP2G ) )
IVEC1(1:IGDRY) = INT( ZVEC1(1:IGDRY) )
ZVEC1(1:IGDRY) = ZVEC1(1:IGDRY) - REAL( IVEC1(1:IGDRY) )
!
ZVEC2(1:IGDRY) = MAX( 1.00001, MIN( REAL(NDRYLBDAR)-0.00001, &
- XDRYINTP1R * LOG( ZVEC2(1:IGDRY) ) + XDRYINTP2R ) )
+ XDRYINTP1R * LOG( ZVECLBDAR(1:IGDRY) ) + XDRYINTP2R ) )
IVEC2(1:IGDRY) = INT( ZVEC2(1:IGDRY) )
ZVEC2(1:IGDRY) = ZVEC2(1:IGDRY) - REAL( IVEC2(1:IGDRY) )
!
@@ -261,19 +254,18 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
- XKER_RDRYG(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) &
* (ZVEC1(JJ) - 1.0)
END DO
- ZZW(:) = 0.
+!
DO JJ = 1, IGDRY
- ZZW(I1(JJ)) = ZVEC3(JJ)
+ JL = I1(JJ)
+ ZZW1(JL,4) = MIN( PRRS(JL),XFRDRYG*ZVEC3(JJ) & ! RRDRYG
+ *( ZVECLBDAR(JJ)**(-4) )*( ZVECLBDAG(JJ)**XCXG ) &
+ *( PRHODREF(JL)**(-XCEXVT-1.) ) &
+ *( XLBRDRYG1/( ZVECLBDAG(JJ)**2 ) + &
+ XLBRDRYG2/( ZVECLBDAG(JJ) * ZVECLBDAR(JJ) ) + &
+ XLBRDRYG3/( ZVECLBDAR(JJ)**2) ) )
END DO
-!
- WHERE( GDRY(:) )
- ZZW1(:,4) = MIN( PRRS(:),XFRDRYG*ZZW(:) & ! RRDRYG
- *( PLBDAR(:)**(-4) )*( PLBDAG(:)**XCXG ) &
- *( PRHODREF(:)**(-XCEXVT-1.) ) &
- *( XLBRDRYG1/( PLBDAG(:)**2 ) + &
- XLBRDRYG2/( PLBDAG(:) * PLBDAR(:) ) + &
- XLBRDRYG3/( PLBDAR(:)**2) ) )
- END WHERE
+ DEALLOCATE(ZVECLBDAR)
+ DEALLOCATE(ZVECLBDAG)
DEALLOCATE(IVEC2)
DEALLOCATE(IVEC1)
DEALLOCATE(ZVEC3)
@@ -285,7 +277,6 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!* 6.3 compute the Wet growth case
!
- ZZW(:) = 0.0
PRWETG(:) = 0.0
WHERE( PRGT(:)>XRTMIN(6) )
ZZW1(:,5) = MIN( PRIS(:), &
@@ -310,7 +301,6 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!* 6.4 Select Wet or Dry case
!
- ZZW(:) = 0.0
IF ( KRR == 7 ) THEN
WHERE( PRGT(:)>XRTMIN(6) .AND. PZT(:)<XTT &
.AND. & ! Wet
@@ -344,14 +334,13 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
WHERE( PRGT(:)>XRTMIN(6) .AND. PZT(:)<XTT &
.AND. & ! Wet
PRDRYG(:)>=PRWETG(:) .AND. PRWETG(:)>0.0 ) ! case
- ZZW(:) = PRWETG(:)
PRCS(:) = PRCS(:) - ZZW1(:,1)
PRIS(:) = PRIS(:) - ZZW1(:,5)
PRSS(:) = PRSS(:) - ZZW1(:,6)
- PRGS(:) = PRGS(:) + ZZW(:)
+ PRGS(:) = PRGS(:) + PRWETG(:)
!
- PRRS(:) = PRRS(:) - ZZW(:) + ZZW1(:,5) + ZZW1(:,6) + ZZW1(:,1)
- PTHS(:) = PTHS(:) + (ZZW(:)-ZZW1(:,5)-ZZW1(:,6))*(PLSFACT(:)-PLVFACT(:))
+ PRRS(:) = PRRS(:) - PRWETG(:) + ZZW1(:,5) + ZZW1(:,6) + ZZW1(:,1)
+ PTHS(:) = PTHS(:) + (PRWETG(:)-ZZW1(:,5)-ZZW1(:,6))*(PLSFACT(:)-PLVFACT(:))
! f(L_f*(RCWETG+RRWETG))
END WHERE
END IF
@@ -417,8 +406,7 @@ REAL, DIMENSION(size(PRHODREF),7) :: ZZW1 ! Work arrays
!
!* 6.5 Melting of the graupeln
!
- ZZW(:) = 0.0
- WHERE( (PRGT(:)>XRTMIN(6)) .AND. (PRGS(:)>0.0) .AND. (PZT(:)>XTT) )
+ WHERE( PRGT(:)>XRTMIN(6) .AND. PRGS(:)>0.0 .AND. PZT(:)>XTT )
ZZW(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
ZZW(:) = PKA(:)*(XTT-PZT(:)) + &
( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
diff --git a/src/MNH/rain_ice_fast_rh.f90 b/src/MNH/rain_ice_fast_rh.f90
index d5f7813dddbf85de69c0662300e7e733718c7efa..cedf7ceb49ac1a4c725839c600fd18d0fb74e4ce 100644
--- a/src/MNH/rain_ice_fast_rh.f90
+++ b/src/MNH/rain_ice_fast_rh.f90
@@ -7,6 +7,7 @@
! P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
! P. Wautelet 03/06/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+! P. Wautelet 05/06/2019: optimisations
!-----------------------------------------------------------------
MODULE MODE_RAIN_ICE_FAST_RH
@@ -72,56 +73,53 @@ REAL, DIMENSION(:), intent(inout) :: PUSW ! Undersaturation over wat
!* 0.2 declaration of local variables
!
INTEGER :: IHAIL, IGWET
-INTEGER :: JJ
-INTEGER, DIMENSION(size(PRHODREF)) :: I1
+INTEGER :: JJ, JL
+INTEGER, DIMENSION(size(PRHODREF)) :: I1H, I1W
INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC1, IVEC2 ! Vectors of indices for interpolations
-LOGICAL, DIMENSION(size(PRHODREF)) :: GWET ! Test where to compute wet growth
-LOGICAL, DIMENSION(size(PRHODREF)) :: GHAIL ! Test where to compute hail growth
REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1,ZVEC2,ZVEC3 ! Work vectors for interpolations
+REAL, DIMENSION(:), ALLOCATABLE :: ZVECLBDAG, ZVECLBDAH, ZVECLBDAS
REAL, DIMENSION(size(PRHODREF)) :: ZZW ! Work array
REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!
!-------------------------------------------------------------------------------
!
IHAIL = 0
- DO JJ = 1, SIZE(GHAIL)
+ DO JJ = 1, SIZE(PRHT)
IF ( PRHT(JJ)>XRTMIN(7) ) THEN
IHAIL = IHAIL + 1
- I1(IHAIL) = JJ
- GHAIL(JJ) = .TRUE.
- ELSE
- GHAIL(JJ) = .FALSE.
+ I1H(IHAIL) = JJ
END IF
END DO
!
IF( IHAIL>0 ) THEN
!
!* 7.2 compute the Wet growth of hail
-!
- WHERE ( GHAIL(:) )
- PLBDAH(:) = XLBH*( PRHODREF(:)*MAX( PRHT(:),XRTMIN(7) ) )**XLBEXH
- END WHERE
!
ZZW1(:,:) = 0.0
- WHERE( GHAIL(:) .AND. ((PRCT(:)>XRTMIN(2) .AND. PRCS(:)>0.0)) )
- ZZW(:) = PLBDAH(:)**(XCXH-XDH-2.0) * PRHODREF(:)**(-XCEXVT)
- ZZW1(:,1) = MIN( PRCS(:),XFWETH * PRCT(:) * ZZW(:) ) ! RCWETH
- END WHERE
- WHERE( GHAIL(:) .AND. ((PRIT(:)>XRTMIN(4) .AND. PRIS(:)>0.0)) )
- ZZW(:) = PLBDAH(:)**(XCXH-XDH-2.0) * PRHODREF(:)**(-XCEXVT)
- ZZW1(:,2) = MIN( PRIS(:),XFWETH * PRIT(:) * ZZW(:) ) ! RIWETH
- END WHERE
+!
+ DO JJ = 1, IHAIL
+ JL = I1H(JJ)
+ PLBDAH(JL) = XLBH * ( PRHODREF(JL) * MAX( PRHT(JL), XRTMIN(7) ) )**XLBEXH
+
+ IF ( PRCT(JL)>XRTMIN(2) .AND. PRCS(JL)>0.0 ) THEN
+ ZZW(JL) = PLBDAH(JL)**(XCXH-XDH-2.0) * PRHODREF(JL)**(-XCEXVT)
+ ZZW1(JL,1) = MIN( PRCS(JL),XFWETH * PRCT(JL) * ZZW(JL) ) ! RCWETH
+ END IF
+
+ IF ( PRIT(JL)>XRTMIN(4) .AND. PRIS(JL)>0.0 ) THEN
+ ZZW(JL) = PLBDAH(JL)**(XCXH-XDH-2.0) * PRHODREF(JL)**(-XCEXVT)
+ ZZW1(JL,2) = MIN( PRIS(JL),XFWETH * PRIT(JL) * ZZW(JL) ) ! RIWETH
+ END IF
+ END DO
!
!* 7.2.1 accretion of aggregates on the hailstones
!
IGWET = 0
- DO JJ = 1, SIZE(GWET)
- IF ( GHAIL(JJ) .AND. PRST(JJ)>XRTMIN(5) .AND. PRSS(JJ)>0.0 ) THEN
+ DO JJ = 1, IHAIL
+ JL = I1H(JJ)
+ IF ( PRST(JL)>XRTMIN(5) .AND. PRSS(JL)>0.0 ) THEN
IGWET = IGWET + 1
- I1(IGWET) = JJ
- GWET(JJ) = .TRUE.
- ELSE
- GWET(JJ) = .FALSE.
+ I1W(IGWET) = JL
END IF
END DO
!
@@ -129,6 +127,8 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!
!* 7.2.2 allocations
!
+ ALLOCATE(ZVECLBDAH(IGWET))
+ ALLOCATE(ZVECLBDAS(IGWET))
ALLOCATE(ZVEC1(IGWET))
ALLOCATE(ZVEC2(IGWET))
ALLOCATE(ZVEC3(IGWET))
@@ -137,22 +137,20 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!
!* 7.2.3 select the (PLBDAH,PLBDAS) couplet
!
- DO JJ = 1, IGWET
- ZVEC1(JJ) = PLBDAH(I1(JJ))
- ZVEC2(JJ) = PLBDAS(I1(JJ))
- END DO
+ ZVECLBDAH(1:IGWET) = PLBDAH(I1W(1:IGWET))
+ ZVECLBDAS(1:IGWET) = PLBDAS(I1W(1:IGWET))
!
!* 7.2.4 find the next lower indice for the PLBDAG and for the PLBDAS
! in the geometrical set of (Lbda_h,Lbda_s) couplet use to
! tabulate the SWETH-kernel
!
ZVEC1(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAH)-0.00001, &
- XWETINTP1H * LOG( ZVEC1(1:IGWET) ) + XWETINTP2H ) )
+ XWETINTP1H * LOG( ZVECLBDAH(1:IGWET) ) + XWETINTP2H ) )
IVEC1(1:IGWET) = INT( ZVEC1(1:IGWET) )
ZVEC1(1:IGWET) = ZVEC1(1:IGWET) - REAL( IVEC1(1:IGWET) )
!
ZVEC2(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAS)-0.00001, &
- XWETINTP1S * LOG( ZVEC2(1:IGWET) ) + XWETINTP2S ) )
+ XWETINTP1S * LOG( ZVECLBDAS(1:IGWET) ) + XWETINTP2S ) )
IVEC2(1:IGWET) = INT( ZVEC2(1:IGWET) )
ZVEC2(1:IGWET) = ZVEC2(1:IGWET) - REAL( IVEC2(1:IGWET) )
!
@@ -167,19 +165,18 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
- XKER_SWETH(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) &
* (ZVEC1(JJ) - 1.0)
END DO
- ZZW(:) = 0.
+!
DO JJ = 1, IGWET
- ZZW(I1(JJ)) = ZVEC3(JJ)
+ JL = I1W(JJ)
+ ZZW1(JL,3) = MIN( PRSS(JL),XFSWETH*ZVEC3(JJ) & ! RSWETH
+ *( ZVECLBDAS(JJ)**(XCXS-XBS) )*( ZVECLBDAH(JJ)**XCXH ) &
+ *( PRHODREF(JL)**(-XCEXVT-1.) ) &
+ *( XLBSWETH1/( ZVECLBDAH(JJ)**2 ) + &
+ XLBSWETH2/( ZVECLBDAH(JJ) * ZVECLBDAS(JJ) ) + &
+ XLBSWETH3/( ZVECLBDAS(JJ)**2) ) )
END DO
-!
- WHERE( GWET(:) )
- ZZW1(:,3) = MIN( PRSS(:),XFSWETH*ZZW(:) & ! RSWETH
- *( PLBDAS(:)**(XCXS-XBS) )*( PLBDAH(:)**XCXH ) &
- *( PRHODREF(:)**(-XCEXVT-1.) ) &
- *( XLBSWETH1/( PLBDAH(:)**2 ) + &
- XLBSWETH2/( PLBDAH(:) * PLBDAS(:) ) + &
- XLBSWETH3/( PLBDAS(:)**2) ) )
- END WHERE
+ DEALLOCATE(ZVECLBDAS)
+ DEALLOCATE(ZVECLBDAH)
DEALLOCATE(IVEC2)
DEALLOCATE(IVEC1)
DEALLOCATE(ZVEC3)
@@ -190,13 +187,11 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!* 7.2.6 accretion of graupeln on the hailstones
!
IGWET = 0
- DO JJ = 1, SIZE(GWET)
- IF ( GHAIL(JJ) .AND. PRGT(JJ)>XRTMIN(6) .AND. PRGS(JJ)>0.0 ) THEN
+ DO JJ = 1, IHAIL
+ JL = I1H(JJ)
+ IF ( PRGT(JL)>XRTMIN(6) .AND. PRGS(JL)>0.0 ) THEN
IGWET = IGWET + 1
- I1(IGWET) = JJ
- GWET(JJ) = .TRUE.
- ELSE
- GWET(JJ) = .FALSE.
+ I1W(IGWET) = JL
END IF
END DO
!
@@ -204,6 +199,8 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!
!* 7.2.7 allocations
!
+ ALLOCATE(ZVECLBDAG(IGWET))
+ ALLOCATE(ZVECLBDAH(IGWET))
ALLOCATE(ZVEC1(IGWET))
ALLOCATE(ZVEC2(IGWET))
ALLOCATE(ZVEC3(IGWET))
@@ -212,22 +209,20 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!
!* 7.2.8 select the (PLBDAH,PLBDAG) couplet
!
- DO JJ = 1, IGWET
- ZVEC1(JJ) = PLBDAH(I1(JJ))
- ZVEC2(JJ) = PLBDAG(I1(JJ))
- END DO
+ ZVECLBDAG(1:IGWET) = PLBDAG(I1W(1:IGWET))
+ ZVECLBDAH(1:IGWET) = PLBDAH(I1W(1:IGWET))
!
!* 7.2.9 find the next lower indice for the PLBDAH and for the PLBDAG
! in the geometrical set of (Lbda_h,Lbda_g) couplet use to
! tabulate the GWETH-kernel
!
ZVEC1(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAG)-0.00001, &
- XWETINTP1H * LOG( ZVEC1(1:IGWET) ) + XWETINTP2H ) )
+ XWETINTP1H * LOG( ZVECLBDAH(1:IGWET) ) + XWETINTP2H ) )
IVEC1(1:IGWET) = INT( ZVEC1(1:IGWET) )
ZVEC1(1:IGWET) = ZVEC1(1:IGWET) - REAL( IVEC1(1:IGWET) )
!
ZVEC2(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAG)-0.00001, &
- XWETINTP1G * LOG( ZVEC2(1:IGWET) ) + XWETINTP2G ) )
+ XWETINTP1G * LOG( ZVECLBDAG(1:IGWET) ) + XWETINTP2G ) )
IVEC2(1:IGWET) = INT( ZVEC2(1:IGWET) )
ZVEC2(1:IGWET) = ZVEC2(1:IGWET) - REAL( IVEC2(1:IGWET) )
!
@@ -242,19 +237,18 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
- XKER_GWETH(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) &
* (ZVEC1(JJ) - 1.0)
END DO
- ZZW(:) = 0.
+!
DO JJ = 1, IGWET
- ZZW(I1(JJ)) = ZVEC3(JJ)
+ JL = I1W(JJ)
+ ZZW1(JL,5) = MAX(MIN( PRGS(JL),XFGWETH*ZVEC3(JJ) & ! RGWETH
+ *( ZVECLBDAG(JJ)**(XCXG-XBG) )*( ZVECLBDAH(JJ)**XCXH ) &
+ *( PRHODREF(JL)**(-XCEXVT-1.) ) &
+ *( XLBGWETH1/( ZVECLBDAH(JJ)**2 ) + &
+ XLBGWETH2/( ZVECLBDAH(JJ) * ZVECLBDAG(JJ) ) + &
+ XLBGWETH3/( ZVECLBDAG(JJ)**2) ) ),0. )
END DO
-!
- WHERE( GWET(:) )
- ZZW1(:,5) = MAX(MIN( PRGS(:),XFGWETH*ZZW(:) & ! RGWETH
- *( PLBDAG(:)**(XCXG-XBG) )*( PLBDAH(:)**XCXH ) &
- *( PRHODREF(:)**(-XCEXVT-1.) ) &
- *( XLBGWETH1/( PLBDAH(:)**2 ) + &
- XLBGWETH2/( PLBDAH(:) * PLBDAG(:) ) + &
- XLBGWETH3/( PLBDAG(:)**2) ) ),0. )
- END WHERE
+ DEALLOCATE(ZVECLBDAH)
+ DEALLOCATE(ZVECLBDAG)
DEALLOCATE(IVEC2)
DEALLOCATE(IVEC1)
DEALLOCATE(ZVEC3)
@@ -264,45 +258,47 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!
!* 7.3 compute the Wet growth of hail
!
- ZZW(:) = 0.0
- WHERE( GHAIL(:) .AND. PZT(:)<XTT )
- ZZW(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
- ZZW(:) = PKA(:)*(XTT-PZT(:)) + &
- ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
- *(XESTT-ZZW(:))/(XRV*PZT(:)) )
+ DO JJ = 1, IHAIL
+ JL = I1H(JJ)
+ IF ( PZT(JL)<XTT ) THEN
+ ZZW(JL) = PRVT(JL)*PPRES(JL)/((XMV/XMD)+PRVT(JL)) ! Vapor pressure
+ ZZW(JL) = PKA(JL)*(XTT-PZT(JL)) + &
+ ( PDV(JL)*(XLVTT + ( XCPV - XCL ) * ( PZT(JL) - XTT )) &
+ *(XESTT-ZZW(JL))/(XRV*PZT(JL)) )
!
! compute RWETH
!
- ZZW(:) = MAX(0., ( ZZW(:) * ( X0DEPH* PLBDAH(:)**XEX0DEPH + &
- X1DEPH*PCJ(:)*PLBDAH(:)**XEX1DEPH ) + &
- ( ZZW1(:,2)+ZZW1(:,3)+ZZW1(:,5) ) * &
- ( PRHODREF(:)*(XLMTT+(XCI-XCL)*(XTT-PZT(:))) ) ) / &
- ( PRHODREF(:)*(XLMTT-XCL*(XTT-PZT(:))) ) )
+ ZZW(JL) = MAX(0., ( ZZW(JL) * ( X0DEPH* PLBDAH(JL)**XEX0DEPH + &
+ X1DEPH*PCJ(JL)*PLBDAH(JL)**XEX1DEPH ) + &
+ ( ZZW1(JL,2)+ZZW1(JL,3)+ZZW1(JL,5) ) * &
+ ( PRHODREF(JL)*(XLMTT+(XCI-XCL)*(XTT-PZT(JL))) ) ) / &
+ ( PRHODREF(JL)*(XLMTT-XCL*(XTT-PZT(JL))) ) )
!
- ZZW1(:,6) = MAX( ZZW(:) - ZZW1(:,2) - ZZW1(:,3) - ZZW1(:,5),0.) ! RCWETH+RRWETH
- END WHERE
- WHERE ( GHAIL(:) .AND. PZT(:)<XTT .AND. ZZW1(:,6)/=0.)
+ ZZW1(JL,6) = MAX( ZZW(JL) - ZZW1(JL,2) - ZZW1(JL,3) - ZZW1(JL,5),0.) ! RCWETH+RRWETH
+ IF ( ZZW1(JL,6)/=0.) THEN
!
! limitation of the available rainwater mixing ratio (RRWETH < RRS !)
!
- ZZW1(:,4) = MAX( 0.0,MIN( ZZW1(:,6),PRRS(:)+ZZW1(:,1) ) )
- PUSW(:) = ZZW1(:,4) / ZZW1(:,6)
- ZZW1(:,2) = ZZW1(:,2)*PUSW(:)
- ZZW1(:,3) = ZZW1(:,3)*PUSW(:)
- ZZW1(:,5) = ZZW1(:,5)*PUSW(:)
- ZZW(:) = ZZW1(:,4) + ZZW1(:,2) + ZZW1(:,3) + ZZW1(:,5)
+ ZZW1(JL,4) = MAX( 0.0,MIN( ZZW1(JL,6),PRRS(JL)+ZZW1(JL,1) ) )
+ PUSW(JL) = ZZW1(JL,4) / ZZW1(JL,6)
+ ZZW1(JL,2) = ZZW1(JL,2)*PUSW(JL)
+ ZZW1(JL,3) = ZZW1(JL,3)*PUSW(JL)
+ ZZW1(JL,5) = ZZW1(JL,5)*PUSW(JL)
+ ZZW(JL) = ZZW1(JL,4) + ZZW1(JL,2) + ZZW1(JL,3) + ZZW1(JL,5)
!
!* 7.1.6 integrate the Wet growth of hail
!
- PRCS(:) = PRCS(:) - ZZW1(:,1)
- PRIS(:) = PRIS(:) - ZZW1(:,2)
- PRSS(:) = PRSS(:) - ZZW1(:,3)
- PRGS(:) = PRGS(:) - ZZW1(:,5)
- PRHS(:) = PRHS(:) + ZZW(:)
- PRRS(:) = MAX( 0.0,PRRS(:) - ZZW1(:,4) + ZZW1(:,1) )
- PTHS(:) = PTHS(:) + ZZW1(:,4)*(PLSFACT(:)-PLVFACT(:))
+ PRCS(JL) = PRCS(JL) - ZZW1(JL,1)
+ PRIS(JL) = PRIS(JL) - ZZW1(JL,2)
+ PRSS(JL) = PRSS(JL) - ZZW1(JL,3)
+ PRGS(JL) = PRGS(JL) - ZZW1(JL,5)
+ PRHS(JL) = PRHS(JL) + ZZW(JL)
+ PRRS(JL) = MAX( 0.0,PRRS(JL) - ZZW1(JL,4) + ZZW1(JL,1) )
+ PTHS(JL) = PTHS(JL) + ZZW1(JL,4)*(PLSFACT(JL)-PLVFACT(JL))
! f(L_f*(RCWETH+RRWETH))
- END WHERE
+ END IF
+ END IF
+ END DO
END IF
IF (LBUDGET_TH) CALL BUDGET ( &
UNPACK(PTHS(:),MASK=OMICRO(:,:,:),FIELD=PTHS3D)*PRHODJ3D(:,:,:),&
@@ -358,24 +354,25 @@ REAL, DIMENSION(size(PRHODREF),6) :: ZZW1 ! Work arrays
!
!* 7.5 Melting of the hailstones
!
- ZZW(:) = 0.0
- WHERE( GHAIL(:) .AND. (PRHS(:)>0.0) .AND. (PZT(:)>XTT) )
- ZZW(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
- ZZW(:) = PKA(:)*(XTT-PZT(:)) + &
- ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
- *(XESTT-ZZW(:))/(XRV*PZT(:)) )
+ DO JJ = 1, IHAIL
+ JL = I1H(JJ)
+ IF( PRHS(JL)>0.0 .AND. PZT(JL)>XTT ) THEN
+ ZZW(JL) = PRVT(JL)*PPRES(JL)/((XMV/XMD)+PRVT(JL)) ! Vapor pressure
+ ZZW(JL) = PKA(JL)*(XTT-PZT(JL)) + &
+ ( PDV(JL)*(XLVTT + ( XCPV - XCL ) * ( PZT(JL) - XTT )) &
+ *(XESTT-ZZW(JL))/(XRV*PZT(JL)) )
!
! compute RHMLTR
!
- ZZW(:) = MIN( PRHS(:), MAX( 0.0,( -ZZW(:) * &
- ( X0DEPH* PLBDAH(:)**XEX0DEPH + &
- X1DEPH*PCJ(:)*PLBDAH(:)**XEX1DEPH ) - &
- ZZW1(:,6)*( PRHODREF(:)*XCL*(XTT-PZT(:))) ) / &
- ( PRHODREF(:)*XLMTT ) ) )
- PRRS(:) = PRRS(:) + ZZW(:)
- PRHS(:) = PRHS(:) - ZZW(:)
- PTHS(:) = PTHS(:) - ZZW(:)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(-RHMLTR))
- END WHERE
+ ZZW(JL) = MIN( PRHS(JL), MAX( 0.0,( -ZZW(JL) * &
+ ( X0DEPH* PLBDAH(JL)**XEX0DEPH + &
+ X1DEPH*PCJ(JL)*PLBDAH(JL)**XEX1DEPH ) ) / &
+ ( PRHODREF(JL)*XLMTT ) ) )
+ PRRS(JL) = PRRS(JL) + ZZW(JL)
+ PRHS(JL) = PRHS(JL) - ZZW(JL)
+ PTHS(JL) = PTHS(JL) - ZZW(JL)*(PLSFACT(JL)-PLVFACT(JL)) ! f(L_f*(-RHMLTR))
+ END IF
+ END DO
END IF
IF (LBUDGET_TH) CALL BUDGET ( &
diff --git a/src/MNH/rain_ice_fast_ri.f90 b/src/MNH/rain_ice_fast_ri.f90
index 2ed51c7ffa059e697affd4a4bce0a4bebb3f32a6..782f79c9eaffbbb699c155fbb85e0872d27e23c7 100644
--- a/src/MNH/rain_ice_fast_ri.f90
+++ b/src/MNH/rain_ice_fast_ri.f90
@@ -5,6 +5,7 @@
!-----------------------------------------------------------------
! Modifications:
! P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
+! P. Wautelet 05/06/2019: optimisations
!-----------------------------------------------------------------
MODULE MODE_RAIN_ICE_FAST_RI
@@ -57,9 +58,7 @@ REAL, DIMENSION(size(PRHODREF)) :: ZZW ! Work array
!
!* 7.1 cloud ice melting
!
- ZZW(:) = 0.0
- WHERE( (PRIS(:)>0.0) .AND. (PZT(:)>XTT) )
- ZZW(:) = PRIS(:)
+ WHERE( PRIS(:)>0.0 .AND. PZT(:)>XTT )
PRCS(:) = PRCS(:) + PRIS(:)
PTHS(:) = PTHS(:) - PRIS(:)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(-RIMLTC))
PRIS(:) = 0.0
@@ -77,9 +76,7 @@ REAL, DIMENSION(size(PRHODREF)) :: ZZW ! Work array
!
!* 7.2 Bergeron-Findeisen effect: RCBERI
!
- ZZW(:) = 0.0
- WHERE( (PRCS(:)>0.0) .AND. (PSSI(:)>0.0) .AND. &
- (PRIT(:)>XRTMIN(4)) .AND. (PCIT(:)>0.0) )
+ WHERE( PRCS(:)>0.0 .AND. PSSI(:)>0.0 .AND. PRIT(:)>XRTMIN(4) .AND. PCIT(:)>0.0 )
ZZW(:) = MIN(1.E8,XLBI*( PRHODREF(:)*PRIT(:)/PCIT(:) )**XLBEXI) ! Lbda_i
ZZW(:) = MIN( PRCS(:),( PSSI(:) / (PRHODREF(:)*PAI(:)) ) * PCIT(:) * &
( X0DEPI/ZZW(:) + X2DEPI*PCJ(:)*PCJ(:)/ZZW(:)**(XDI+2.0) ) )
diff --git a/src/MNH/rain_ice_fast_rs.f90 b/src/MNH/rain_ice_fast_rs.f90
index 7f90da8c02c112e140808145d38b1024f21cc37c..5f5f9713eba6e28af53732ca458c168ad098f5fb 100644
--- a/src/MNH/rain_ice_fast_rs.f90
+++ b/src/MNH/rain_ice_fast_rs.f90
@@ -7,6 +7,7 @@
! P. Wautelet 25/02/2019: split rain_ice (cleaner and easier to maintain/debug)
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
! P. Wautelet 03/06/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
+! P. Wautelet 05/06/2019: optimisations
!-----------------------------------------------------------------
MODULE MODE_RAIN_ICE_FAST_RS
@@ -69,28 +70,22 @@ REAL, DIMENSION(:), INTENT(INOUT) :: PTHS ! Theta source
!* 0.2 declaration of local variables
!
INTEGER :: IGRIM, IGACC
-INTEGER :: JJ
+INTEGER :: JJ, JL
INTEGER, DIMENSION(size(PRHODREF)) :: I1
INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC1, IVEC2 ! Vectors of indices for interpolations
-LOGICAL, DIMENSION(size(PRHODREF)) :: GRIM ! Test where to compute riming
-LOGICAL, DIMENSION(size(PRHODREF)) :: GACC ! Test where to compute accretion
REAL, DIMENSION(size(PRHODREF)) :: ZZW ! Work array
REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1,ZVEC2,ZVEC3 ! Work vectors for interpolations
-REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
+REAL, DIMENSION(:), ALLOCATABLE :: ZVECLBDAR, ZVECLBDAS
+REAL, DIMENSION(:), ALLOCATABLE :: ZZW1, ZZW2, ZZW3, ZZW4 ! Work arrays
!-------------------------------------------------------------------------------
!
!* 5.1 cloud droplet riming of the aggregates
-!
- ZZW1(:,:) = 0.0
!
IGRIM = 0
- DO JJ = 1, SIZE(GRIM)
- IF (PRCT(JJ)>XRTMIN(2) .AND. PRST(JJ)>XRTMIN(5) .AND. PRCS(JJ)>0.0 .AND. PZT(JJ)<XTT ) THEN
+ DO JJ = 1, SIZE(PRCT)
+ IF ( PRCT(JJ)>XRTMIN(2) .AND. PRST(JJ)>XRTMIN(5) .AND. PRCS(JJ)>0.0 .AND. PZT(JJ)<XTT ) THEN
IGRIM = IGRIM + 1
I1(IGRIM) = JJ
- GRIM(JJ) = .TRUE.
- ELSE
- GRIM(JJ) = .FALSE.
END IF
END DO
!
@@ -98,22 +93,24 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
!
! 5.1.0 allocations
!
+ ALLOCATE(ZVECLBDAS(IGRIM))
ALLOCATE(ZVEC1(IGRIM))
ALLOCATE(ZVEC2(IGRIM))
ALLOCATE(IVEC2(IGRIM))
+ ALLOCATE(ZZW1(IGRIM))
+ ALLOCATE(ZZW2(IGRIM))
+ ALLOCATE(ZZW3(IGRIM))
!
! 5.1.1 select the PLBDAS
!
- DO JJ = 1, IGRIM
- ZVEC1(JJ) = PLBDAS(I1(JJ))
- END DO
+ ZVECLBDAS(1:IGRIM) = PLBDAS(I1(1:IGRIM))
!
! 5.1.2 find the next lower indice for the PLBDAS in the geometrical
! set of Lbda_s used to tabulate some moments of the incomplete
! gamma function
!
ZVEC2(1:IGRIM) = MAX( 1.00001, MIN( REAL(NGAMINC)-0.00001, &
- XRIMINTP1 * LOG( ZVEC1(1:IGRIM) ) + XRIMINTP2 ) )
+ XRIMINTP1 * LOG( ZVECLBDAS(1:IGRIM) ) + XRIMINTP2 ) )
IVEC2(1:IGRIM) = INT( ZVEC2(1:IGRIM) )
ZVEC2(1:IGRIM) = ZVEC2(1:IGRIM) - REAL( IVEC2(1:IGRIM) )
!
@@ -122,53 +119,53 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
!
ZVEC1(1:IGRIM) = XGAMINC_RIM1( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM) &
- XGAMINC_RIM1( IVEC2(1:IGRIM) )*(ZVEC2(1:IGRIM) - 1.0)
- ZZW(:) = 0.
- DO JJ = 1, IGRIM
- ZZW(I1(JJ)) = ZVEC1(JJ)
- END DO
!
! 5.1.4 riming of the small sized aggregates
!
- WHERE ( GRIM(:) )
- ZZW1(:,1) = MIN( PRCS(:), &
- XCRIMSS * ZZW(:) * PRCT(:) & ! RCRIMSS
- * PLBDAS(:)**XEXCRIMSS &
- * PRHODREF(:)**(-XCEXVT) )
- PRCS(:) = PRCS(:) - ZZW1(:,1)
- PRSS(:) = PRSS(:) + ZZW1(:,1)
- PTHS(:) = PTHS(:) + ZZW1(:,1)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(RCRIMSS))
- END WHERE
+ DO JJ = 1, IGRIM
+ JL = I1(JJ)
+ ZZW1(JJ) = MIN( PRCS(JL), &
+ XCRIMSS * ZVEC1(JJ) * PRCT(JL) & ! RCRIMSS
+ * ZVECLBDAS(JJ)**XEXCRIMSS &
+ * PRHODREF(JL)**(-XCEXVT) )
+ PRCS(JL) = PRCS(JL) - ZZW1(JJ)
+ PRSS(JL) = PRSS(JL) + ZZW1(JJ)
+ PTHS(JL) = PTHS(JL) + ZZW1(JJ)*(PLSFACT(JL)-PLVFACT(JL)) ! f(L_f*(RCRIMSS))
+ END DO
!
! 5.1.5 perform the linear interpolation of the normalized
! "XBS"-moment of the incomplete gamma function
!
ZVEC1(1:IGRIM) = XGAMINC_RIM2( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM) &
- XGAMINC_RIM2( IVEC2(1:IGRIM) )*(ZVEC2(1:IGRIM) - 1.0)
- ZZW(:) = 0.
- DO JJ = 1, IGRIM
- ZZW(I1(JJ)) = ZVEC1(JJ)
- END DO
!
! 5.1.6 riming-conversion of the large sized aggregates into graupeln
!
!
- WHERE ( GRIM(:) .AND. (PRSS(:)>0.0) )
- ZZW1(:,2) = MIN( PRCS(:), &
- XCRIMSG * PRCT(:) & ! RCRIMSG
- * PLBDAS(:)**XEXCRIMSG &
- * PRHODREF(:)**(-XCEXVT) &
- - ZZW1(:,1) )
- ZZW1(:,3) = MIN( PRSS(:), &
- XSRIMCG * PLBDAS(:)**XEXSRIMCG & ! RSRIMCG
- * (1.0 - ZZW(:) )/(PTSTEP*PRHODREF(:)) )
- PRCS(:) = PRCS(:) - ZZW1(:,2)
- PRSS(:) = PRSS(:) - ZZW1(:,3)
- PRGS(:) = PRGS(:) + ZZW1(:,2)+ZZW1(:,3)
- PTHS(:) = PTHS(:) + ZZW1(:,2)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(RCRIMSG))
- END WHERE
+ DO JJ = 1, IGRIM
+ JL = I1(JJ)
+ IF ( PRSS(JL) > 0.0 ) THEN
+ ZZW2(JJ) = MIN( PRCS(JL), &
+ XCRIMSG * PRCT(JL) & ! RCRIMSG
+ * ZVECLBDAS(JJ)**XEXCRIMSG &
+ * PRHODREF(JL)**(-XCEXVT) &
+ - ZZW1(JJ) )
+ ZZW3(JJ) = MIN( PRSS(JL), &
+ XSRIMCG * ZVECLBDAS(JJ)**XEXSRIMCG & ! RSRIMCG
+ * (1.0 - ZVEC1(JJ) )/(PTSTEP*PRHODREF(JL)) )
+ PRCS(JL) = PRCS(JL) - ZZW2(JJ)
+ PRSS(JL) = PRSS(JL) - ZZW3(JJ)
+ PRGS(JL) = PRGS(JL) + ZZW2(JJ)+ZZW3(JJ)
+ PTHS(JL) = PTHS(JL) + ZZW2(JJ)*(PLSFACT(JL)-PLVFACT(JL)) ! f(L_f*(RCRIMSG))
+ END IF
+ END DO
+ DEALLOCATE(ZZW3)
+ DEALLOCATE(ZZW2)
+ DEALLOCATE(ZZW1)
DEALLOCATE(IVEC2)
DEALLOCATE(ZVEC2)
DEALLOCATE(ZVEC1)
+ DEALLOCATE(ZVECLBDAS)
END IF
IF (LBUDGET_TH) CALL BUDGET ( &
UNPACK(PTHS(:),MASK=OMICRO(:,:,:),FIELD=PTHS3D)*PRHODJ3D(:,:,:), &
@@ -185,16 +182,11 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
!
!* 5.2 rain accretion onto the aggregates
!
- ZZW1(:,2:3) = 0.0
-
IGACC = 0
- DO JJ = 1, SIZE(GACC)
- IF (PRRT(JJ)>XRTMIN(3) .AND. PRST(JJ)>XRTMIN(5) .AND. PRRS(JJ)>0.0 .AND. PZT(JJ)<XTT ) THEN
+ DO JJ = 1, SIZE(PRRT)
+ IF ( PRRT(JJ)>XRTMIN(3) .AND. PRST(JJ)>XRTMIN(5) .AND. PRRS(JJ)>0.0 .AND. PZT(JJ)<XTT ) THEN
IGACC = IGACC + 1
I1(IGACC) = JJ
- GACC(JJ) = .TRUE.
- ELSE
- GACC(JJ) = .FALSE.
END IF
END DO
!
@@ -202,30 +194,33 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
!
! 5.2.0 allocations
!
+ ALLOCATE(ZVECLBDAR(IGACC))
+ ALLOCATE(ZVECLBDAS(IGACC))
ALLOCATE(ZVEC1(IGACC))
ALLOCATE(ZVEC2(IGACC))
ALLOCATE(ZVEC3(IGACC))
ALLOCATE(IVEC1(IGACC))
ALLOCATE(IVEC2(IGACC))
+ ALLOCATE(ZZW2(IGACC))
+ ALLOCATE(ZZW3(IGACC))
+ ALLOCATE(ZZW4(IGACC))
!
! 5.2.1 select the (PLBDAS,PLBDAR) couplet
!
- DO JJ = 1, IGACC
- ZVEC1(JJ) = PLBDAS(I1(JJ))
- ZVEC2(JJ) = PLBDAR(I1(JJ))
- END DO
+ ZVECLBDAS(1:IGACC) = PLBDAS(I1(1:IGACC))
+ ZVECLBDAR(1:IGACC) = PLBDAR(I1(1:IGACC))
!
! 5.2.2 find the next lower indice for the PLBDAS and for the PLBDAR
! in the geometrical set of (Lbda_s,Lbda_r) couplet use to
! tabulate the RACCSS-kernel
!
ZVEC1(1:IGACC) = MAX( 1.00001, MIN( REAL(NACCLBDAS)-0.00001, &
- XACCINTP1S * LOG( ZVEC1(1:IGACC) ) + XACCINTP2S ) )
+ XACCINTP1S * LOG( ZVECLBDAS(1:IGACC) ) + XACCINTP2S ) )
IVEC1(1:IGACC) = INT( ZVEC1(1:IGACC) )
ZVEC1(1:IGACC) = ZVEC1(1:IGACC) - REAL( IVEC1(1:IGACC) )
!
ZVEC2(1:IGACC) = MAX( 1.00001, MIN( REAL(NACCLBDAR)-0.00001, &
- XACCINTP1R * LOG( ZVEC2(1:IGACC) ) + XACCINTP2R ) )
+ XACCINTP1R * LOG( ZVECLBDAR(1:IGACC) ) + XACCINTP2R ) )
IVEC2(1:IGACC) = INT( ZVEC2(1:IGACC) )
ZVEC2(1:IGACC) = ZVEC2(1:IGACC) - REAL( IVEC2(1:IGACC) )
!
@@ -240,24 +235,21 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
- XKER_RACCSS(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) &
* (ZVEC1(JJ) - 1.0)
END DO
- ZZW(:) = 0.
- DO JJ = 1, IGACC
- ZZW(I1(JJ)) = ZVEC3(JJ)
- END DO
!
! 5.2.4 raindrop accretion on the small sized aggregates
!
- WHERE ( GACC(:) )
- ZZW1(:,2) = & !! coef of RRACCS
- XFRACCSS*( PLBDAS(:)**XCXS )*( PRHODREF(:)**(-XCEXVT-1.) ) &
- *( XLBRACCS1/((PLBDAS(:)**2) ) + &
- XLBRACCS2/( PLBDAS(:) * PLBDAR(:) ) + &
- XLBRACCS3/( (PLBDAR(:)**2)) )/PLBDAR(:)**4
- ZZW1(:,4) = MIN( PRRS(:),ZZW1(:,2)*ZZW(:) ) ! RRACCSS
- PRRS(:) = PRRS(:) - ZZW1(:,4)
- PRSS(:) = PRSS(:) + ZZW1(:,4)
- PTHS(:) = PTHS(:) + ZZW1(:,4)*(PLSFACT(:)-PLVFACT(:)) ! f(L_f*(RRACCSS))
- END WHERE
+ DO JJ = 1, IGACC
+ JL = I1(JJ)
+ ZZW2(JJ) = & !! coef of RRACCS
+ XFRACCSS*( ZVECLBDAS(JJ)**XCXS )*( PRHODREF(JL)**(-XCEXVT-1.) ) &
+ *( XLBRACCS1/((ZVECLBDAS(JJ)**2) ) + &
+ XLBRACCS2/( ZVECLBDAS(JJ) * ZVECLBDAR(JJ) ) + &
+ XLBRACCS3/( (ZVECLBDAR(JJ)**2)) )/ZVECLBDAR(JJ)**4
+ ZZW4(JJ) = MIN( PRRS(JL),ZZW2(JJ)*ZVEC3(JJ) ) ! RRACCSS
+ PRRS(JL) = PRRS(JL) - ZZW4(JJ)
+ PRSS(JL) = PRSS(JL) + ZZW4(JJ)
+ PTHS(JL) = PTHS(JL) + ZZW4(JJ)*(PLSFACT(JL)-PLVFACT(JL)) ! f(L_f*(RRACCSS))
+ END DO
!
! 5.2.4b perform the bilinear interpolation of the normalized
! RACCS-kernel
@@ -271,7 +263,7 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
* (ZVEC2(JJ) - 1.0)
END DO
DO JJ = 1, IGACC
- ZZW1(I1(JJ), 2) = ZZW1(I1(JJ), 2 ) * ZVEC3(JJ)
+ ZZW2(JJ) = ZZW2(JJ) * ZVEC3(JJ)
END DO
!! RRACCS!
! 5.2.5 perform the bilinear interpolation of the normalized
@@ -285,34 +277,38 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
- XKER_SACCRG(IVEC2(JJ) ,IVEC1(JJ) )*(ZVEC1(JJ) - 1.0) ) &
* (ZVEC2(JJ) - 1.0)
END DO
- ZZW(:) = 0.
- DO JJ = 1, IGACC
- ZZW(I1(JJ)) = ZVEC3(JJ)
- END DO
!
! 5.2.6 raindrop accretion-conversion of the large sized aggregates
! into graupeln
!
- WHERE ( GACC(:) .AND. (PRSS(:)>0.0) )
- ZZW1(:,2) = MAX( MIN( PRRS(:),ZZW1(:,2)-ZZW1(:,4) ),0.0 ) ! RRACCSG
- END WHERE
- WHERE ( GACC(:) .AND. (PRSS(:)>0.0) .AND. ZZW1(:,2)>0.0 )
- ZZW1(:,3) = MIN( PRSS(:),XFSACCRG*ZZW(:)* & ! RSACCRG
- ( PLBDAS(:)**(XCXS-XBS) )*( PRHODREF(:)**(-XCEXVT-1.) ) &
- *( XLBSACCR1/((PLBDAR(:)**2) ) + &
- XLBSACCR2/( PLBDAR(:) * PLBDAS(:) ) + &
- XLBSACCR3/( (PLBDAS(:)**2)) )/PLBDAR(:) )
- PRRS(:) = PRRS(:) - ZZW1(:,2)
- PRSS(:) = PRSS(:) - ZZW1(:,3)
- PRGS(:) = PRGS(:) + ZZW1(:,2)+ZZW1(:,3)
- PTHS(:) = PTHS(:) + ZZW1(:,2)*(PLSFACT(:)-PLVFACT(:)) !
- ! f(L_f*(RRACCSG))
- END WHERE
+ DO JJ = 1, IGACC
+ JL = I1(JJ)
+ IF ( PRSS(JL) > 0.0 ) THEN
+ ZZW2(JJ) = MAX( MIN( PRRS(JL),ZZW2(JJ)-ZZW4(JJ) ),0.0 ) ! RRACCSG
+ IF ( ZZW2(JJ) > 0.0 ) THEN
+ ZZW3(JJ) = MIN( PRSS(JL),XFSACCRG*ZVEC3(JJ)* & ! RSACCRG
+ ( ZVECLBDAS(JJ)**(XCXS-XBS) )*( PRHODREF(JL)**(-XCEXVT-1.) ) &
+ *( XLBSACCR1/((ZVECLBDAR(JJ)**2) ) + &
+ XLBSACCR2/( ZVECLBDAR(JJ) * ZVECLBDAS(JJ) ) + &
+ XLBSACCR3/( (ZVECLBDAS(JJ)**2)) )/ZVECLBDAR(JJ) )
+ PRRS(JL) = PRRS(JL) - ZZW2(JJ)
+ PRSS(JL) = PRSS(JL) - ZZW3(JJ)
+ PRGS(JL) = PRGS(JL) + ZZW2(JJ)+ZZW3(JJ)
+ PTHS(JL) = PTHS(JL) + ZZW2(JJ)*(PLSFACT(JL)-PLVFACT(JL)) !
+ ! f(L_f*(RRACCSG))
+ END IF
+ END IF
+ END DO
+ DEALLOCATE(ZZW4)
+ DEALLOCATE(ZZW3)
+ DEALLOCATE(ZZW2)
DEALLOCATE(IVEC2)
DEALLOCATE(IVEC1)
DEALLOCATE(ZVEC3)
DEALLOCATE(ZVEC2)
DEALLOCATE(ZVEC1)
+ DEALLOCATE(ZVECLBDAS)
+ DEALLOCATE(ZVECLBDAR)
END IF
IF (LBUDGET_TH) CALL BUDGET ( &
UNPACK(PTHS(:),MASK=OMICRO(:,:,:),FIELD=PTHS3D)*PRHODJ3D(:,:,:), &
@@ -329,8 +325,7 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
!
!* 5.3 Conversion-Melting of the aggregates
!
- ZZW(:) = 0.0
- WHERE( (PRST(:)>XRTMIN(5)) .AND. (PRSS(:)>0.0) .AND. (PZT(:)>XTT) )
+ WHERE( PRST(:)>XRTMIN(5) .AND. PRSS(:)>0.0 .AND. PZT(:)>XTT )
ZZW(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
ZZW(:) = PKA(:)*(XTT-PZT(:)) + &
( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PZT(:) - XTT )) &
@@ -340,9 +335,7 @@ REAL, DIMENSION(size(PRHODREF),4) :: ZZW1 ! Work arrays
!
ZZW(:) = MIN( PRSS(:), XFSCVMG*MAX( 0.0,( -ZZW(:) * &
( X0DEPS* PLBDAS(:)**XEX0DEPS + &
- X1DEPS*PCJ(:)*PLBDAS(:)**XEX1DEPS ) - &
- ( ZZW1(:,1)+ZZW1(:,4) ) * &
- ( PRHODREF(:)*XCL*(XTT-PZT(:))) ) / &
+ X1DEPS*PCJ(:)*PLBDAS(:)**XEX1DEPS ) ) / &
( PRHODREF(:)*XLMTT ) ) )
!
! note that RSCVMG = RSMLT*XFSCVMG but no heat is exchanged (at the rate RSMLT)