diff --git a/src/ICE/icedyn_adv_pra.F90 b/src/ICE/icedyn_adv_pra.F90
index 27b828ec5aba945d8547c7b213b0959025231cea..763f71a6611c1f296982b755a5aeafb9b3de820e 100644
--- a/src/ICE/icedyn_adv_pra.F90
+++ b/src/ICE/icedyn_adv_pra.F90
@@ -391,158 +391,152 @@ CONTAINS
DO jl = 1, jcat ! loop on categories
!
! Limitation of moments.
- DO jj = Njs0 - jj0, Nje0 + jj0
-
- DO ji = Nis0 - 1, Nie0 + 1
-
- zpsm = psm (ji,jj,jl) ! optimization
- zps0 = ps0 (ji,jj,jl)
- zpsx = psx (ji,jj,jl)
- zpsxx = psxx(ji,jj,jl)
- zpsy = psy (ji,jj,jl)
- zpsyy = psyy(ji,jj,jl)
- zpsxy = psxy(ji,jj,jl)
-
- ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
- zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
- !
- zslpmax = MAX( 0._wp, zps0 )
- zs1max = 1.5 * zslpmax
- zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) )
- zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
- rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
-
- zps0 = zslpmax
- zpsx = zs1new * rswitch
- zpsxx = zs2new * rswitch
- zpsy = zpsy * rswitch
- zpsyy = zpsyy * rswitch
- zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
-
- ! Calculate fluxes and moments between boxes i<-->i+1
- ! ! Flux from i to i+1 WHEN u GT 0
- zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
- zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
- zalfq = zalf * zalf
- zalf1 = 1.0 - zalf
- zalf1q = zalf1 * zalf1
- !
- zfm (ji,jj) = zalf * zpsm
- zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
- zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx )
- zfxx(ji,jj) = zalf * zpsxx * zalfq
- zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy )
- zfxy(ji,jj) = zalfq * zpsxy
- zfyy(ji,jj) = zalf * zpsyy
-
- ! ! Readjust moments remaining in the box.
- zpsm = zpsm - zfm(ji,jj)
- zps0 = zps0 - zf0(ji,jj)
- zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
- zpsxx = zalf1 * zalf1q * zpsxx
- zpsy = zpsy - zfy (ji,jj)
- zpsyy = zpsyy - zfyy(ji,jj)
- zpsxy = zalf1q * zpsxy
- !
- psm (ji,jj,jl) = zpsm ! optimization
- ps0 (ji,jj,jl) = zps0
- psx (ji,jj,jl) = zpsx
- psxx(ji,jj,jl) = zpsxx
- psy (ji,jj,jl) = zpsy
- psyy(ji,jj,jl) = zpsyy
- psxy(ji,jj,jl) = zpsxy
- !
- END DO
-
- DO ji = Nis0 - 1, Nie0
- ! ! Flux from i+1 to i when u LT 0.
- zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
- zalg (ji,jj) = zalf
- zalfq = zalf * zalf
- zalf1 = 1.0 - zalf
- zalg1 (ji,jj) = zalf1
- zalf1q = zalf1 * zalf1
- zalg1q(ji,jj) = zalf1q
- !
- zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
- zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
- & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
- zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
- zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq
- zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
- zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
- zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
- END DO
-
- DO ji = Nis0, Nie0
- !
- zpsm = psm (ji,jj,jl) ! optimization
- zps0 = ps0 (ji,jj,jl)
- zpsx = psx (ji,jj,jl)
- zpsxx = psxx(ji,jj,jl)
- zpsy = psy (ji,jj,jl)
- zpsyy = psyy(ji,jj,jl)
- zpsxy = psxy(ji,jj,jl)
- ! ! Readjust moments remaining in the box.
- zbt = zbet(ji-1,jj)
- zbt1 = 1.0 - zbet(ji-1,jj)
- !
- zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
- zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
- zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
- zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
- zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) )
- zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
- zpsxy = zalg1q(ji-1,jj) * zpsxy
+ DO_2D( 1, 1, jj0, jj0 )
+ !
+ zpsm = psm (ji,jj,jl) ! optimization
+ zps0 = ps0 (ji,jj,jl)
+ zpsx = psx (ji,jj,jl)
+ zpsxx = psxx(ji,jj,jl)
+ zpsy = psy (ji,jj,jl)
+ zpsyy = psyy(ji,jj,jl)
+ zpsxy = psxy(ji,jj,jl)
+
+ ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
+ zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
+ !
+ zslpmax = MAX( 0._wp, zps0 )
+ zs1max = 1.5 * zslpmax
+ zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) )
+ zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
+ rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
+
+ zps0 = zslpmax
+ zpsx = zs1new * rswitch
+ zpsxx = zs2new * rswitch
+ zpsy = zpsy * rswitch
+ zpsyy = zpsyy * rswitch
+ zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+
+ ! Calculate fluxes and moments between boxes i<-->i+1
+ ! ! Flux from i to i+1 WHEN u GT 0
+ zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
+ zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
+ zalfq = zalf * zalf
+ zalf1 = 1.0 - zalf
+ zalf1q = zalf1 * zalf1
+ !
+ zfm (ji,jj) = zalf * zpsm
+ zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
+ zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx )
+ zfxx(ji,jj) = zalf * zpsxx * zalfq
+ zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy )
+ zfxy(ji,jj) = zalfq * zpsxy
+ zfyy(ji,jj) = zalf * zpsyy
+
+ ! ! Readjust moments remaining in the box.
+ zpsm = zpsm - zfm(ji,jj)
+ zps0 = zps0 - zf0(ji,jj)
+ zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
+ zpsxx = zalf1 * zalf1q * zpsxx
+ zpsy = zpsy - zfy (ji,jj)
+ zpsyy = zpsyy - zfyy(ji,jj)
+ zpsxy = zalf1q * zpsxy
+ !
+ psm (ji,jj,jl) = zpsm ! optimization
+ ps0 (ji,jj,jl) = zps0
+ psx (ji,jj,jl) = zpsx
+ psxx(ji,jj,jl) = zpsxx
+ psy (ji,jj,jl) = zpsy
+ psyy(ji,jj,jl) = zpsyy
+ psxy(ji,jj,jl) = zpsxy
+ !
+ END_2D
- ! Put the temporary moments into appropriate neighboring boxes.
- ! ! Flux from i to i+1 IF u GT 0.
- zbt = zbet(ji-1,jj)
- zbt1 = 1.0 - zbet(ji-1,jj)
- zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm
- zalf = zbt * zfm(ji-1,jj) / zpsm
- zalf1 = 1.0 - zalf
- ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
- !
- zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0
- zpsx = zbt * ( ( zalf * zfx(ji-1,jj) + zalf1 * zpsx ) + 3.0 * ztemp ) + zbt1 * zpsx
- zpsxx = zbt * ( ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx ) &
- & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
- & + zbt1 * zpsxx
- zpsxy = zbt * ( ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy ) &
- & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * zpsy ) ) &
- & + zbt1 * zpsxy
- zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy
- zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy
+ DO_2D( 1, 0, jj0, jj0 )
+ ! ! Flux from i+1 to i when u LT 0.
+ zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
+ zalg (ji,jj) = zalf
+ zalfq = zalf * zalf
+ zalf1 = 1.0 - zalf
+ zalg1 (ji,jj) = zalf1
+ zalf1q = zalf1 * zalf1
+ zalg1q(ji,jj) = zalf1q
+ !
+ zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
+ zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
+ & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
+ zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
+ zfxx (ji,jj) = zfxx(ji,jj) + zalf * ( psxx(ji+1,jj,jl) * zalfq )
+ zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
+ zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
+ zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
+ END_2D
- ! ! Flux from i+1 to i IF u LT 0.
- zbt = zbet(ji,jj)
- zbt1 = 1.0 - zbet(ji,jj)
- zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
- zalf = zbt1 * zfm(ji,jj) / zpsm
- zalf1 = 1.0 - zalf
- ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
- !
- zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
- zpsx = zbt * zpsx + zbt1 * ( ( zalf * zfx(ji,jj) + zalf1 * zpsx ) + 3.0 * ztemp )
- zpsxx = zbt * zpsxx + zbt1 * ( ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx ) &
- & + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) &
- & + ( zalf1 - zalf ) * ztemp ) )
- zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj) + zalf1 * zpsxy ) &
- & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) )
- zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) )
- zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) )
- !
- psm (ji,jj,jl) = zpsm ! optimization
- ps0 (ji,jj,jl) = zps0
- psx (ji,jj,jl) = zpsx
- psxx(ji,jj,jl) = zpsxx
- psy (ji,jj,jl) = zpsy
- psyy(ji,jj,jl) = zpsyy
- psxy(ji,jj,jl) = zpsxy
- END DO
+ DO_2D( 0, 0, jj0, jj0 )
!
- END DO
+ zpsm = psm (ji,jj,jl) ! optimization
+ zps0 = ps0 (ji,jj,jl)
+ zpsx = psx (ji,jj,jl)
+ zpsxx = psxx(ji,jj,jl)
+ zpsy = psy (ji,jj,jl)
+ zpsyy = psyy(ji,jj,jl)
+ zpsxy = psxy(ji,jj,jl)
+ ! ! Readjust moments remaining in the box.
+ zbt = zbet(ji-1,jj)
+ zbt1 = 1.0 - zbet(ji-1,jj)
+ !
+ zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
+ zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
+ zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
+ zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
+ zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) )
+ zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
+ zpsxy = zalg1q(ji-1,jj) * zpsxy
+
+ ! Put the temporary moments into appropriate neighboring boxes.
+ ! ! Flux from i to i+1 IF u GT 0.
+ zbt = zbet(ji-1,jj)
+ zbt1 = 1.0 - zbet(ji-1,jj)
+ zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji-1,jj) )
+ zalf = zbt * zfm(ji-1,jj) / zpsm
+ zalf1 = 1.0 - zalf
+ ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
+ !
+ zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji-1,jj) )
+ zpsx = zbt1 * zpsx + zbt * ( ( zalf * zfx(ji-1,jj) + zalf1 * zpsx ) + 3.0 * ztemp )
+ zpsxx = zbt1 * zpsxx + zbt * ( ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx ) &
+ & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy ) &
+ & + 3.0 * ( -zalf1 * zfy(ji-1,jj) + zalf * zpsy ) )
+ zpsy = zbt1 * zpsy + zbt * ( zpsy + zfy (ji-1,jj) )
+ zpsyy = zbt1 * zpsyy + zbt * ( zpsyy + zfyy(ji-1,jj) )
+
+ ! ! Flux from i+1 to i IF u LT 0.
+ zbt = zbet(ji,jj)
+ zbt1 = 1.0 - zbet(ji,jj)
+ zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji ,jj) )
+ zalf = zbt1 * zfm(ji ,jj) / zpsm
+ zalf1 = 1.0 - zalf
+ ztemp = -zalf * zps0 + zalf1 * zf0(ji ,jj)
+ !
+ zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji ,jj) )
+ zpsx = zbt * zpsx + zbt1 * ( ( zalf * zfx(ji ,jj) + zalf1 * zpsx ) + 3.0 * ztemp )
+ zpsxx = zbt * zpsxx + zbt1 * ( ( zalf * zalf * zfxx(ji ,jj) + zalf1 * zalf1 * zpsxx ) &
+ & + 5.0 * ( zalf * zalf1 * ( -zpsx + zfx(ji ,jj) ) + ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji ,jj) + zalf1 * zpsxy ) &
+ & + 3.0 * ( zalf1 * zfy(ji ,jj) - zalf * zpsy ) )
+ zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji ,jj) )
+ zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji ,jj) )
+ !
+ psm (ji,jj,jl) = zpsm ! optimization
+ ps0 (ji,jj,jl) = zps0
+ psx (ji,jj,jl) = zpsx
+ psxx(ji,jj,jl) = zpsxx
+ psy (ji,jj,jl) = zpsy
+ psyy(ji,jj,jl) = zpsyy
+ psxy(ji,jj,jl) = zpsxy
+ !
+ END_2D
!
END DO
!
@@ -602,14 +596,14 @@ CONTAINS
zslpmax = MAX( 0._wp, zps0 )
zs1max = 1.5 * zslpmax
zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) )
- zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) )
+ zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new ) - zslpmax, zpsyy ) )
rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
!
zps0 = zslpmax
- zpsx = zpsx * rswitch
- zpsxx = zpsxx * rswitch
- zpsy = zs1new * rswitch
- zpsyy = zs2new * rswitch
+ zpsx = zpsx * rswitch
+ zpsxx = zpsxx * rswitch
+ zpsy = zs1new * rswitch
+ zpsyy = zs2new * rswitch
zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
! Calculate fluxes and moments between boxes j<-->j+1
@@ -620,19 +614,19 @@ CONTAINS
zalf1 = 1.0 - zalf
zalf1q = zalf1 * zalf1
!
- zfm (ji,jj) = zalf * zpsm
- zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) )
- zfy (ji,jj) = zalfq *( zpsy + 3.0*zalf1*zpsyy )
- zfyy(ji,jj) = zalf * zalfq * zpsyy
- zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy )
- zfxy(ji,jj) = zalfq * zpsxy
- zfxx(ji,jj) = zalf * zpsxx
+ zfm (ji,jj) = zalf * zpsm
+ zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1 - zalf) * zpsyy ) )
+ zfy (ji,jj) = zalfq * ( zpsy + 3.0 * zalf1 * zpsyy )
+ zfyy(ji,jj) = zalf * zpsyy * zalfq
+ zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy )
+ zfxy(ji,jj) = zalfq * zpsxy
+ zfxx(ji,jj) = zalf * zpsxx
!
! ! Readjust moments remaining in the box.
zpsm = zpsm - zfm(ji,jj)
zps0 = zps0 - zf0(ji,jj)
zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy )
- zpsyy = zalf1 * zalf1q * zpsyy
+ zpsyy = zalf1 * zalf1q * zpsyy
zpsx = zpsx - zfx(ji,jj)
zpsxx = zpsxx - zfxx(ji,jj)
zpsxy = zalf1q * zpsxy
@@ -660,16 +654,13 @@ CONTAINS
zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) &
& - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
- zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq
+ zfyy (ji,jj) = zfyy(ji,jj) + zalf * ( psyy(ji,jj+1,jl) * zalfq )
zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl)
zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl)
END_2D
DO_2D( ji0, ji0, 0, 0 )
- ! ! Readjust moments remaining in the box.
- zbt = zbet(ji,jj-1)
- zbt1 = ( 1.0 - zbet(ji,jj-1) )
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
@@ -678,53 +669,52 @@ CONTAINS
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
+ ! ! Readjust moments remaining in the box.
+ zbt = zbet(ji,jj-1)
+ zbt1 = ( 1.0 - zbet(ji,jj-1) )
!
- zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
- zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
+ zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
+ zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
- zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) )
- zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
+ zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) )
+ zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
zpsxy = zalg1q(ji,jj-1) * zpsxy
! Put the temporary moments into appropriate neighboring boxes.
! ! Flux from j to j+1 IF v GT 0.
zbt = zbet(ji,jj-1)
zbt1 = 1.0 - zbet(ji,jj-1)
- zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm
- zalf = zbt * zfm(ji,jj-1) / zpsm
+ zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji,jj-1) )
+ zalf = zbt * zfm(ji,jj-1) / zpsm
zalf1 = 1.0 - zalf
- ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
+ ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
!
- zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0
- zpsy = zbt * ( ( zalf * zfy(ji,jj-1) + zalf1 * zpsy ) + 3.0 * ztemp ) &
- & + zbt1 * zpsy
- zpsyy = zbt * ( ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy ) &
- & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &
- & + zbt1 * zpsyy
- zpsxy = zbt * ( ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy ) &
- & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) &
- & + zbt1 * zpsxy
- zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx
- zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx
+ zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji,jj-1) )
+ zpsy = zbt1 * zpsy + zbt * ( ( zalf * zfy(ji,jj-1) + zalf1 * zpsy ) + 3.0 * ztemp )
+ zpsyy = zbt1 * zpsyy + zbt * ( ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy ) &
+ & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy ) &
+ & + 3.0 * ( -zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )
+ zpsx = zbt1 * zpsx + zbt * ( zpsx + zfx (ji,jj-1) )
+ zpsxx = zbt1 * zpsxx + zbt * ( zpsxx + zfxx(ji,jj-1) )
! ! Flux from j+1 to j IF v LT 0.
zbt = zbet(ji,jj)
zbt1 = 1.0 - zbet(ji,jj)
- zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
- zalf = zbt1 * zfm(ji,jj) / zpsm
+ zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
+ zalf = zbt1 * zfm(ji,jj) / zpsm
zalf1 = 1.0 - zalf
- ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
+ ztemp = -zalf * zps0 + zalf1 * zf0(ji,jj)
!
- zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
- zpsy = zbt * zpsy + zbt1 * ( ( zalf * zfy(ji,jj) + zalf1 * zpsy ) + 3.0 * ztemp )
- zpsyy = zbt * zpsyy + zbt1 * ( ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy ) &
- & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) &
- & + ( zalf1 - zalf ) * ztemp ) )
- zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj) + zalf1 * zpsxy ) &
- & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) )
- zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) )
- zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) )
+ zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj ) )
+ zpsy = zbt * zpsy + zbt1 * ( ( zalf * zfy(ji,jj ) + zalf1 * zpsy ) + 3.0 * ztemp )
+ zpsyy = zbt * zpsyy + zbt1 * ( ( zalf * zalf * zfyy(ji,jj ) + zalf1 * zalf1 * zpsyy ) &
+ & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj ) ) + ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj ) + zalf1 * zpsxy ) &
+ & + 3.0 * ( zalf1 * zfx(ji,jj ) - zalf * zpsx ) )
+ zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj ) )
+ zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj ) )
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0