! Generated by TAPENADE (INRIA, Tropics team)
! Version 2.1 - (Id: 1.33 vmp Stable - Wed Jan 12 13:35:07 MET 2005)
!
! $Id: ice_transport_mpdata.F,v 1.11 2004/02/25 17:35:14 eclare Exp $
!=======================================================================
!BOP
!
! !MODULE: ice_transport_mpdata - horizontal advection (via mpdata)
!
! !DESCRIPTION:
!
! Calculates horizontal advection using mpdata (Multidimensional
! Positive Definite Advection Transport Algorithm).
!
! !REVISION HISTORY:
!
! author Elizabeth C. Hunke
!
! !INTERFACE:
!
MODULE ICE_TRANSPORT_MPDATA
USE ice_calendar
USE ice_constants
USE ice_domain
USE ice_exit
USE ice_fileunits
USE ice_flux
USE ice_grid
USE ice_itd
USE ice_kinds_mod
USE ice_model_size
USE ice_state
USE ice_timers
USE ice_transport_remap
USE ice_work
CONTAINS
SUBROUTINE TRANSPORT_MPDATA()
IMPLICIT NONE
INCLUDE 'DIFFSIZES.inc'
! Hint: nbdirsmax should be the maximum number of differentiation directions
INTEGER(KIND=INT_KIND),PARAMETER :: narr=1+5*ncat
INTEGER(KIND=INT_KIND) :: i, j, k, n, narrays
REAL(KIND=DBL_KIND) :: worke(imt_local, jmt_local, ntilay)
REAL(KIND=DBL_KIND) :: works(imt_local, jmt_local, narr)
!
! !DESCRIPTION:
!
! Computes the transport equations for one timestep using mpdata. Sets
! several fields into a work array and passes it to mpdata routine.
!
! !REVISION HISTORY:
!
! author Elizabeth C. Hunke
!
! !USES:
!
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! number of state variable arrays
!
! standard indices
! counter for number of state variable arrays
!
!
CALL ICE_TIMER_START(3)
! advection
!-----------------------------------------------------------------
! Make sure state variables are not too close to zero.
!-----------------------------------------------------------------
!
!-----------------------------------------------------------------
! fill work arrays with fields to be advected
!-----------------------------------------------------------------
! two arrays are used for performance (balance memory/cache vs
! number of bound calls); one array or more than two may perform
! better depending on the machine used, number of processors, etc.
! --tested on SGI R2000, using 4 pes for the ice model under MPI
!-----------------------------------------------------------------
CALL CHECK_STATE()
!
!
DO j=jlo,jhi
DO i=ilo,ihi
works(i, j, 1) = aice0(i, j)
END DO
END DO
!
narrays = 1
DO n=1,ncat
DO j=jlo,jhi
DO i=ilo,ihi
works(i, j, narrays+1) = aicen(i, j, n)
works(i, j, narrays+2) = vicen(i, j, n)
works(i, j, narrays+3) = vsnon(i, j, n)
works(i, j, narrays+4) = -(aicen(i, j, n)*tsfcn(i, j, n))
! for Tsnow<0
works(i, j, narrays+5) = -(esnon(i, j, n)+rhos*lfresh*vsnon(i&
& , j, n))
END DO
END DO
! if there were 3 arrays in loop, use 3 instead and change narr
!
narrays = narrays + 5
END DO
!
DO k=1,ntilay
DO j=jlo,jhi
DO i=ilo,ihi
worke(i, j, k) = -eicen(i, j, k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------
! update local domain boundaries
!-----------------------------------------------------------------
IF (narr .NE. narrays) WRITE(nu_diag, *) &
& 'Wrong number of arrays in transport bound call'
!
!
CALL BOUND_NARR(narr, works)
CALL BOUND_NARR(ntilay, worke)
!-----------------------------------------------------------------
! advect
!-----------------------------------------------------------------
!
!
CALL MPDATA(narr, works)
CALL MPDATA(ntilay, worke)
!-----------------------------------------------------------------
! retrieve advected fields from work array
!-----------------------------------------------------------------
!
!
DO j=1,jmt_local
DO i=1,imt_local
aice0(i, j) = works(i, j, 1)
END DO
END DO
!
! aice0 is first array
narrays = 1
DO n=1,ncat
DO j=1,jmt_local
DO i=1,imt_local
aicen(i, j, n) = works(i, j, narrays+1)
vicen(i, j, n) = works(i, j, narrays+2)
vsnon(i, j, n) = works(i, j, narrays+3)
IF (aicen(i, j, n) .GT. puny) THEN
tsfcn(i, j, n) = -(works(i, j, narrays+4)/aicen(i, j, n))
ELSE
tsfcn(i, j, n) = tf(i, j)
END IF
! for Tsnow<0
esnon(i, j, n) = -works(i, j, narrays+5) - rhos*lfresh*vsnon(i&
& , j, n)
END DO
END DO
!
narrays = narrays + 5
END DO
!
DO k=1,ntilay
DO j=1,jmt_local
DO i=1,imt_local
eicen(i, j, k) = -worke(i, j, k)
END DO
END DO
END DO
!
CALL ICE_TIMER_STOP(3)
! advection
!-----------------------------------------------------------------
! mask
!-----------------------------------------------------------------
!
DO n=1,ncat
DO j=1,jmt_local
DO i=1,imt_local
IF (.NOT.tmask(i, j)) THEN
aicen(i, j, n) = c0
vicen(i, j, n) = c0
vsnon(i, j, n) = c0
tsfcn(i, j, n) = c0
esnon(i, j, n) = c0
END IF
END DO
END DO
END DO
!
DO k=1,ntilay
DO j=1,jmt_local
DO i=1,imt_local
IF (.NOT.tmask(i, j)) eicen(i, j, k) = c0
END DO
END DO
END DO
END SUBROUTINE TRANSPORT_MPDATA
SUBROUTINE MPDATA(narrays, phi)
IMPLICIT NONE
!
! !DESCRIPTION:
!
! Smolarkiewicz, P. K., 1984: A fully multidimensional positive
! definite advection transport algorithm with small implicit
! diffusion, J. Comput. Phys., 54, 325-362.
!
! !REVISION HISTORY:
!
! author Elizabeth C. Hunke
! Spring 2003: upwind algorithm modified by William Lipscomb (LANL)
! to improve performance.
!
! !USES:
!
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!
!
!EOP
!
! 2nd order MPDATA paramter
!
! function
!
!
!
INCLUDE 'DIFFSIZES.inc'
! Hint: nbdirsmax should be the maximum number of differentiation directions
INTEGER(KIND=INT_KIND),INTENT(IN) :: narrays
REAL(KIND=DBL_KIND), DIMENSION(imt_local, jmt_local, narrays)&
& ,INTENT(INOUT) :: phi
INTEGER(KIND=INT_KIND),PARAMETER :: iord=3
REAL(KIND=DBL_KIND) :: a, dive(ilo:ihi, jlo:jhi), divn(ilo:ihi, jlo:&
& jhi), eps, h, phiavg(imt_local, jmt_local), uee(imt_local, &
& jmt_local, narrays), upwind, vnn(imt_local, jmt_local, narrays), y1&
& , y2
REAL :: abs1
REAL :: abs10
REAL :: abs11
REAL :: abs12
REAL :: abs13
REAL :: abs14
REAL :: abs2
REAL :: abs3
REAL :: abs4
REAL :: abs5
REAL :: abs6
REAL :: abs7
REAL :: abs8
REAL :: abs9
INTEGER(KIND=INT_KIND) :: ikim, jkim
INTEGER(KIND=INT_KIND) :: i, ix, iy, j, k, n
REAL(KIND=DBL_KIND) :: uv_kim(imt_local, jmt_local)
INTRINSIC ABS
!
eps = eps15
! upwind or mpdata
!
IF (advection .EQ. 'upwind') THEN
!
jkim = jlo - 1
DO j=1,jmt_local
jkim = jkim + 1
ikim = ilo - 1
DO i=1,imt_local
ikim = ikim + 1
uv_kim(i, j) = p5*(uvel(ikim, jkim)+uvel(ikim, jkim-1))
END DO
END DO
CALL BOUND(uv_kim)
jkim = jlo - 1
DO j=1,jmt_local
jkim = jkim + 1
ikim = ilo - 1
DO i=1,imt_local
ikim = ikim + 1
uee(ikim, jkim, 1) = uv_kim(i, j)
END DO
END DO
!
jkim = jlo - 1
DO j=1,jmt_local
jkim = jkim + 1
ikim = ilo - 1
DO i=1,imt_local
ikim = ikim + 1
uv_kim(i, j) = p5*(vvel(ikim, jkim)+vvel(ikim-1, jkim))
END DO
END DO
CALL BOUND(uv_kim)
jkim = jlo - 1
DO j=1,jmt_local
jkim = jkim + 1
ikim = ilo - 1
DO i=1,imt_local
ikim = ikim + 1
vnn(ikim, jkim, 1) = uv_kim(i, j)
END DO
END DO
!
!jkim do j=jlo,jhi
!jkim do i=ilo,ihi
!jkim uee(i,j,1)=p5*(uvel(i,j)+uvel(i,j-1))
!jkim vnn(i,j,1)=p5*(vvel(i,j)+vvel(i-1,j))
!jkim enddo
!jkim enddo
!
!jkim call bound(uee(:,:,1))
!jkim call bound(vnn(:,:,1))
!
DO n=1,narrays
DO j=1,jhi
DO i=1,ihi
IF (uee(i, j, 1) .GE. 0.) THEN
abs1 = uee(i, j, 1)
ELSE
abs1 = -uee(i, j, 1)
END IF
IF (uee(i, j, 1) .GE. 0.) THEN
abs9 = uee(i, j, 1)
ELSE
abs9 = -uee(i, j, 1)
END IF
work_l1(i, j) = p5*dt*hte(i, j)*((uee(i, j, 1)+abs1)*phi(i, &
& j, n)+(uee(i, j, 1)-abs9)*phi(i+1, j, n))
IF (vnn(i, j, 1) .GE. 0.) THEN
abs2 = vnn(i, j, 1)
ELSE
abs2 = -vnn(i, j, 1)
END IF
IF (vnn(i, j, 1) .GE. 0.) THEN
abs10 = vnn(i, j, 1)
ELSE
abs10 = -vnn(i, j, 1)
END IF
work_l2(i, j) = p5*dt*htn(i, j)*((vnn(i, j, 1)+abs2)*phi(i, &
& j, n)+(vnn(i, j, 1)-abs10)*phi(i, j+1, n))
END DO
END DO
DO j=jlo,jhi
DO i=ilo,ihi
phi(i, j, n) = phi(i, j, n) - (work_l1(i, j)-work_l1(i-1, j)&
& +work_l2(i, j)-work_l2(i, j-1))/tarea(i, j)
END DO
END DO
END DO
ELSE
! narrays
!
! mpdata with antidiffusive corrections
!
DO n=1,narrays
DO j=jlo,jhi
DO i=ilo,ihi
uee(i, j, n) = p5*(uvel(i, j)+uvel(i, j-1))
vnn(i, j, n) = p5*(vvel(i, j)+vvel(i-1, j))
END DO
END DO
END DO
! narrays
!
CALL BOUND_NARR(narrays, uee)
CALL BOUND_NARR(narrays, vnn)
! upwind
!
DO n=1,narrays
DO j=1,jhi
DO i=1,ihi
IF (uee(i, j, n) .GE. 0.) THEN
abs3 = uee(i, j, n)
ELSE
abs3 = -uee(i, j, n)
END IF
IF (uee(i, j, n) .GE. 0.) THEN
abs11 = uee(i, j, n)
ELSE
abs11 = -uee(i, j, n)
END IF
work_l1(i, j) = p5*dt*hte(i, j)*((uee(i, j, n)+abs3)*phi(i, &
& j, n)+(uee(i, j, n)-abs11)*phi(i+1, j, n))
IF (vnn(i, j, n) .GE. 0.) THEN
abs4 = vnn(i, j, n)
ELSE
abs4 = -vnn(i, j, n)
END IF
IF (vnn(i, j, n) .GE. 0.) THEN
abs12 = vnn(i, j, n)
ELSE
abs12 = -vnn(i, j, n)
END IF
work_l2(i, j) = p5*dt*htn(i, j)*((vnn(i, j, n)+abs4)*phi(i, &
& j, n)+(vnn(i, j, n)-abs12)*phi(i, j+1, n))
END DO
END DO
DO j=jlo,jhi
DO i=ilo,ihi
phi(i, j, n) = phi(i, j, n) - (work_l1(i, j)-work_l1(i-1, j)&
& +work_l2(i, j)-work_l2(i, j-1))/tarea(i, j)
END DO
END DO
END DO
! narrays
!
CALL BOUND_NARR(narrays, phi)
!
! 2nd order MPDATA
DO k=1,iord
!
DO n=1,narrays
!
DO j=1,jhi
DO i=1,ihi
phiavg(i, j) = p25*(phi(i, j, n)+phi(i+1, j, n)+phi(i+1, j&
& +1, n)+phi(i, j+1, n))
END DO
END DO
!
DO j=jlo,jhi
DO i=ilo,ihi
dive(i, j) = ((dyt(i+1, j)*(uee(i+1, j, n)+uee(i, j, n))*&
& phi(i+1, j, n)-dyt(i, j)*(uee(i-1, j, n)+uee(i, j, n))*&
& phi(i, j, n))/(phi(i+1, j, n)+phi(i, j, n)+eps)+(dxu(i&
& , j)*(vnn(i+1, j, n)+vnn(i, j, n))*phiavg(i, j)-dxu(i, &
& j-1)*(vnn(i+1, j-1, n)+vnn(i, j-1, n))*phiavg(i, j-1))/&
& (phiavg(i, j)+phiavg(i, j-1)+eps))/(hte(i, j)*(dxu(i, j&
& )+dxu(i, j-1)))
!
divn(i, j) = ((dxt(i, j+1)*(vnn(i, j+1, n)+vnn(i, j, n))*&
& phi(i, j+1, n)-dxt(i, j)*(vnn(i, j-1, n)+vnn(i, j, n))*&
& phi(i, j, n))/(phi(i, j+1, n)+phi(i, j, n)+eps)+(dyu(i&
& , j)*(uee(i, j+1, n)+uee(i, j, n))*phiavg(i, j)-dyu(i-1&
& , j)*(uee(i-1, j, n)+uee(i-1, j+1, n))*phiavg(i-1, j))/&
& (phiavg(i, j)+phiavg(i-1, j)+eps))/(htn(i, j)*(dyu(i, j&
& )+dyu(i-1, j)))
END DO
END DO
! antidiffusive velocities
!
DO j=jlo,jhi
DO i=ilo,ihi
IF (uee(i, j, n) .GE. 0.) THEN
abs5 = uee(i, j, n)
ELSE
abs5 = -uee(i, j, n)
END IF
uee(i, j, n) = abs5*(phi(i+1, j, n)-phi(i, j, n))/(phi(i+1&
& , j, n)+phi(i, j, n)+eps) - dt*uee(i, j, n)*dive(i, j)
IF (vnn(i, j, n) .GE. 0.) THEN
abs6 = vnn(i, j, n)
ELSE
abs6 = -vnn(i, j, n)
END IF
!
vnn(i, j, n) = abs6*(phi(i, j+1, n)-phi(i, j, n))/(phi(i, &
& j+1, n)+phi(i, j, n)+eps) - dt*vnn(i, j, n)*divn(i, j)
END DO
END DO
END DO
! narrays
!
CALL BOUND_NARR(narrays, uee)
CALL BOUND_NARR(narrays, vnn)
! upwind with antidiffusive velocities
!
DO n=1,narrays
!
DO j=1,jhi
DO i=1,ihi
IF (uee(i, j, n) .GE. 0.) THEN
abs7 = uee(i, j, n)
ELSE
abs7 = -uee(i, j, n)
END IF
IF (uee(i, j, n) .GE. 0.) THEN
abs13 = uee(i, j, n)
ELSE
abs13 = -uee(i, j, n)
END IF
work_l1(i, j) = p5*dt*hte(i, j)*((uee(i, j, n)+abs7)*phi(i&
& , j, n)+(uee(i, j, n)-abs13)*phi(i+1, j, n))
IF (vnn(i, j, n) .GE. 0.) THEN
abs8 = vnn(i, j, n)
ELSE
abs8 = -vnn(i, j, n)
END IF
IF (vnn(i, j, n) .GE. 0.) THEN
abs14 = vnn(i, j, n)
ELSE
abs14 = -vnn(i, j, n)
END IF
work_l2(i, j) = p5*dt*htn(i, j)*((vnn(i, j, n)+abs8)*phi(i&
& , j, n)+(vnn(i, j, n)-abs14)*phi(i, j+1, n))
END DO
END DO
!
ix = -1
DO j=jlo,jhi
DO i=ilo,ihi
phi(i, j, n) = phi(i, j, n) - (work_l1(i, j)-work_l1(i-1, &
& j)+work_l2(i, j)-work_l2(i, j-1))/tarea(i, j)
!
IF (phi(i, j, n) .LT. -eps12) THEN
ix = i
iy = j
ELSE IF (phi(i, j, n) .LT. 0.) THEN
phi(i, j, n) = c0
END IF
END DO
END DO
!
IF (ix .GE. 0) THEN
WRITE(nu_diag, *) my_task, ix, iy, ' transport unstable ', &
& istep, k
WRITE(nu_diag, *) 'mpdata phi = ', 
' n = ', &
& n
CALL ABORT_ICE('(mpdata) transport unstable')
END IF
END DO
!
! narrays
!
IF (k .LT. iord) CALL BOUND_NARR(narrays, phi)
END DO
END IF
END SUBROUTINE MPDATA
END MODULE ICE_TRANSPORT_MPDATA