Skip to content
Snippets Groups Projects
Commit 11dbcd7f authored by WAUTELET Philippe's avatar WAUTELET Philippe
Browse files

Philippe 04/06/2019: optimisation on rain_ice_fast_rg/h/i/s

parent ac24b04b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/MNH/rain_ice_fast_rg.f90
+ 44
− 56
View file @ 11dbcd7f
......@@ -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 )) &
......
src/MNH/rain_ice_fast_rh.f90
+ 107
− 110
View file @ 11dbcd7f
......@@ -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 ( &
......
src/MNH/rain_ice_fast_ri.f90
+ 3
− 6
View file @ 11dbcd7f
......@@ -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) ) )
......
src/MNH/rain_ice_fast_rs.f90
+ 92
− 99
View file @ 11dbcd7f
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme