compare01/src/inc_tracer.F

C $Id: inc_tracer.F,v 1.3 1998/06/15 05:13:54 cnh Exp $
#include "CPP_OPTIONS.h"
#include "CPP_MACROS.h"
C     Procedure name: INC_TRACER
C           Function: Controls steppping forward tracer field.
C            Written: 10th January 1994, by Chris Hill M.I.T.
C            Updated: Fall 1995, Moto Nakamura, M.I.T., mixing schemes
C            Updated: Fall 1995, Curtis Heisey, M.I.T., generalized routine
C           Comments:
CStartofinterface
      SUBROUTINE INC_TRACER (
     I    U, V, W,
     U    Q, QT, QTNM1,
     I    Kredigm,K31,K32,K33,
     I    surfProf,tracerStepping,tracerDiffusion,tracerBhDiffusion,
     I    tracerAdvection,tracerForcing,a4TracerXY,a2TracerXY,
     I    a2TracerZ 
     O    )
      IMPLICIT NONE
C      ======== Global data ============================
#include "SIZE.h"
#include "OPERATORS.h"
#include "GRID.h"
#include "PARAMS.h"
#include "MASKS.h"
C     ======== Routine arguments ======================
      REAL U          (_I3(Nz,Nx,Ny))
      REAL V          (_I3(Nz,Nx,Ny))
      REAL W          (_I3(Nz,Nx,Ny))
      REAL Q  (_I3(Nz,Nx,Ny))
      REAL QT    (_I3(Nz,Nx,Ny))
      REAL QTNM1 (_I3(Nz,Nx,Ny))
      REAL Kredigm (_I3(Nz,Nx,Ny))
      REAL K31 (_I3(Nz,Nx,Ny))
      REAL K32 (_I3(Nz,Nx,Ny))
      REAL K33 (_I3(Nz,Nx,Ny))
      REAL surfProf(Nx,Ny)
      LOGICAL tracerStepping
      LOGICAL tracerDiffusion
      LOGICAL tracerBhDiffusion
      LOGICAL tracerAdvection
      LOGICAL tracerForcing
      REAL a4TracerXY
      REAL a2TracerXY
      REAL a2TracerZ
CEndofinterface
C     Local variables
C     qXM1 - Q shifted -1 in X.
C     qYM1 - Q shifted -1 in Y.
      REAL qXM1      (_I3(Nz,Nx,Ny))
      REAL qYM1      (_I3(Nz,Nx,Ny))
      REAL qZM1      (_I3(Nz,Nx,Ny))
      REAL fluxWest  (_I3(Nz,Nx,Ny))
      REAL fluxSouth (_I3(Nz,Nx,Ny))
      REAL fluxUpper (_I3(Nz,Nx,Ny))
      REAL dqDx (_I3(Nz,Nx,Ny))
      REAL dqDy (_I3(Nz,Nx,Ny))
      REAL dqDz (_I3(Nz,Nx,Ny))
      REAL dqDzBz(_I3(Nz,Nx,Ny))
      integer K

      qXM1 = Cshift(Q,DIM=_I3X,SHIFT=-1)
      qYM1 = Cshift(Q,DIM=_I3Y,SHIFT=-1)
      qZM1 = Cshift(Q,DIM=_I3Z,SHIFT=-1)
      qT        = 0.

      fluxWest  = 0.
      fluxSouth = 0.
      fluxUpper = 0.
      IF ( tracerDiffusion ) THEN
C       Include del**2 3d dissipation of Q, -1.*div[k2.grad(Q)],
C       and/or del**4 lateral dissipation of Q, divh[k4.gradh(divh[gradh(Q))].
        dQdx = (Q-qXM1)*pDdxOp
        dQdy = (Q-qYM1)*pDdyOp
        do K=2,Nz
         dQdz(_I3(K,:,:))=(G*RONIL)*rDZatW(K)*(qZM1(_I3(K,:,:))-Q(_I3(K,:,:)))
        enddo
! Extrapolate downward from above at W=0 surfaces
        dQdz=WMASK*dQdz+(1.-WMASK)*CSHIFT(dQdz,DIM=_I3Z,SHIFT=-1)
! Extrapolate upward from below at Z=0
        dQdz(_I3(1,:,:))=dQdz(_I3(2,:,:))
! Interpolate to tracer points
        dQdzBz=0.5*(dqdz+CSHIFT(dQdz,DIM=_I3Z,SHIFT=+1))*PMASK
        dQdzBz(_I3(Nz,:,:))=dQdz(_I3(Nz-1,:,:))

!??     fluxWest = XA*(Q-qXM1)*pDdxOp*umask
!??     fluxSouth = YA*(Q-qYM1)*pDdyOp*vmask
!??     IF ( tracerBhDiffusion ) fluxUpper=rPVol*PMASK*(
!??  &  Cshift(fluxWest,DIM=_I3X,SHIFT=+1)
!??  &         -fluxWest+Cshift(fluxSouth,DIM=_I3Y,SHIFT=+1)-fluxSouth)
        IF ( tracerBhDiffusion ) THEN
         fluxWest = XA*(Q-qXM1)*pDdxOp*umask
         fluxSouth = YA*(Q-qYM1)*pDdyOp*vmask
         fluxUpper=rPVol*PMASK*( Cshift(fluxWest,DIM=_I3X,SHIFT=+1)
     &         -fluxWest+Cshift(fluxSouth,DIM=_I3Y,SHIFT=+1)-fluxSouth)
        ENDIF

        fluxWest = -UMASK*XA*(
     &    (0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3X,SHIFT=-1))+a2TracerXY)
     &            *dQDx
     &                       )

        fluxSouth = -VMASK*YA*(
     &    (0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3Y,SHIFT=-1))+a2TracerXY)
     &            *dQDy
     &                       )

        IF ( tracerBhDiffusion ) THEN
          fluxWest = fluxWest  +
     &    a4TracerXY*XA*(fluxUpper-Cshift(fluxUpper,DIM=_I3X,SHIFT=-1))*pDdxOp
          fluxSouth = fluxSouth+
     &    a4TracerXY*YA*(fluxUpper-Cshift(fluxUpper,DIM=_I3Y,SHIFT=-1))*pDdyOp
        ENDIF

! Interpolate d/dx Q and d/dy Q to W points
        dQdx=0.5d0*(dQdx+CSHIFT(dQdx,DIM=_I3X,SHIFT=+1))*PMASK  ! Interp to T
        dQdx=0.5d0*(dQdx+CSHIFT(dQdx,DIM=_I3Z,SHIFT=-1))*WMASK  ! Interp to W
        dQdy=0.5d0*(dQdy+CSHIFT(dQdy,DIM=_I3Y,SHIFT=+1))*PMASK  ! Interp to T
        dQdy=0.5d0*(dQdy+CSHIFT(dQdy,DIM=_I3Z,SHIFT=-1))*WMASK  ! Interp to W

!!      fluxUpper = -ZA*0.5*(Kredi+CSHIFT(Kredi,DIM=_I3Z,SHIFT=-1))
!!   &        *( K31*dQdx + K32*dQdy + K33*dQdz )

        fluxUpper = -ZA*(
     &         0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3Z,SHIFT=-1))
     &                  *( K31*dQdx + K32*dQdy )
     &        +0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3Z,SHIFT=-1))
     &                  *( K33*dQdz )   
     &     )

! Add in vertical diffusion (Av d^2/dz^2)
! fluxUpper=-wmask*g*ronil*fluxUpper-a2TracerZ*wmask*ZA*(Q-qZM1))*pDdzOP
        fluxUpper = -WMASK*( (G*RONIL)*fluxUpper-a2TracerZ*ZA*(qZM1-Q)*pDdzOP )
        fluxUpper(_I3(1,:,:)) = 0.
      ENDIF
      IF ( tracerAdvection ) THEN
C      Include transport of temperature by flow field.
       fluxWest  = fluxWest +XA*(qXM1*0.5+Q*0.5)*U*UMASK
       fluxSouth = fluxSouth+YA*(qYM1*0.5+Q*0.5)*V*VMASK
       fluxUpper = fluxUpper+ZA*(qZM1*0.5+Q*0.5)*W*WMASK
      ENDIF
      IF ( tracerDiffusion .OR. tracerAdvection ) THEN
        qT = qT - (
     &              +Cshift(fluxWest,DIM=_I3X,SHIFT=+1)-fluxWest
     &              +Cshift(fluxSouth,DIM=_I3Y,SHIFT=+1)-fluxSouth
     &              -fluxUpper+Cshift(fluxUpper,DIM=_I3Z,SHIFT=+1)
     &            )*rPVol
      ENDIF
      qT(_I3(1,:,:)) = qT(_I3(1,:,:))
     &               + freeSurfFac*W(_I3(1,:,:))*Q(_I3(1,:,:))/delps(1)

      IF(tracerForcing)THEN
C     Include external forcing ( e.g. surface heat fluxes, boundary fluxes )
C      Switch in the appropriate tracer forcing module 
       qT(_I3(1,:,:)) = qT(_I3(1,:,:)) 
     & - 1./(30.*oneDay)*(Q(_I3(1,:,:))-surfProf(:,:))
c    & - surfProf(_I3(2,:,:))
      ENDIF
      qT = qT*PMASK

C
      Q = Q + ( (1.5+abEps)*qT-(0.5+abEps)*qTNM1 )*delT*asyncfac
C

      RETURN
      END
1	C $Id: inc_tracer.F,v 1.3 1998/06/15 05:13:54 cnh Exp $
2	#include "CPP_OPTIONS.h"
3	#include "CPP_MACROS.h"
4	C Procedure name: INC_TRACER
5	C Function: Controls steppping forward tracer field.
6	C Written: 10th January 1994, by Chris Hill M.I.T.
7	C Updated: Fall 1995, Moto Nakamura, M.I.T., mixing schemes
8	C Updated: Fall 1995, Curtis Heisey, M.I.T., generalized routine
9	C Comments:
10	CStartofinterface
11	SUBROUTINE INC_TRACER (
12	I U, V, W,
13	U Q, QT, QTNM1,
14	I Kredigm,K31,K32,K33,
15	I surfProf,tracerStepping,tracerDiffusion,tracerBhDiffusion,
16	I tracerAdvection,tracerForcing,a4TracerXY,a2TracerXY,
17	I a2TracerZ
18	O )
19	IMPLICIT NONE
20	C ======== Global data ============================
21	#include "SIZE.h"
22	#include "OPERATORS.h"
23	#include "GRID.h"
24	#include "PARAMS.h"
25	#include "MASKS.h"
26	C ======== Routine arguments ======================
27	REAL U (_I3(Nz,Nx,Ny))
28	REAL V (_I3(Nz,Nx,Ny))
29	REAL W (_I3(Nz,Nx,Ny))
30	REAL Q (_I3(Nz,Nx,Ny))
31	REAL QT (_I3(Nz,Nx,Ny))
32	REAL QTNM1 (_I3(Nz,Nx,Ny))
33	REAL Kredigm (_I3(Nz,Nx,Ny))
34	REAL K31 (_I3(Nz,Nx,Ny))
35	REAL K32 (_I3(Nz,Nx,Ny))
36	REAL K33 (_I3(Nz,Nx,Ny))
37	REAL surfProf(Nx,Ny)
38	LOGICAL tracerStepping
39	LOGICAL tracerDiffusion
40	LOGICAL tracerBhDiffusion
41	LOGICAL tracerAdvection
42	LOGICAL tracerForcing
43	REAL a4TracerXY
44	REAL a2TracerXY
45	REAL a2TracerZ
46	CEndofinterface
47	C Local variables
48	C qXM1 - Q shifted -1 in X.
49	C qYM1 - Q shifted -1 in Y.
50	REAL qXM1 (_I3(Nz,Nx,Ny))
51	REAL qYM1 (_I3(Nz,Nx,Ny))
52	REAL qZM1 (_I3(Nz,Nx,Ny))
53	REAL fluxWest (_I3(Nz,Nx,Ny))
54	REAL fluxSouth (_I3(Nz,Nx,Ny))
55	REAL fluxUpper (_I3(Nz,Nx,Ny))
56	REAL dqDx (_I3(Nz,Nx,Ny))
57	REAL dqDy (_I3(Nz,Nx,Ny))
58	REAL dqDz (_I3(Nz,Nx,Ny))
59	REAL dqDzBz(_I3(Nz,Nx,Ny))
60	integer K
61
62	qXM1 = Cshift(Q,DIM=_I3X,SHIFT=-1)
63	qYM1 = Cshift(Q,DIM=_I3Y,SHIFT=-1)
64	qZM1 = Cshift(Q,DIM=_I3Z,SHIFT=-1)
65	qT = 0.
66
67	fluxWest = 0.
68	fluxSouth = 0.
69	fluxUpper = 0.
70	IF ( tracerDiffusion ) THEN
71	C Include del*2 3d dissipation of Q, -1.div[k2.grad(Q)],
72	C and/or del**4 lateral dissipation of Q, divh[k4.gradh(divh[gradh(Q))].
73	dQdx = (Q-qXM1)*pDdxOp
74	dQdy = (Q-qYM1)*pDdyOp
75	do K=2,Nz
76	dQdz(_I3(K,:,:))=(GRONIL)rDZatW(K)*(qZM1(_I3(K,:,:))-Q(_I3(K,:,:)))
77	enddo
78	! Extrapolate downward from above at W=0 surfaces
79	dQdz=WMASKdQdz+(1.-WMASK)CSHIFT(dQdz,DIM=_I3Z,SHIFT=-1)
80	! Extrapolate upward from below at Z=0
81	dQdz(_I3(1,:,:))=dQdz(_I3(2,:,:))
82	! Interpolate to tracer points
83	dQdzBz=0.5(dqdz+CSHIFT(dQdz,DIM=_I3Z,SHIFT=+1))PMASK
84	dQdzBz(_I3(Nz,:,:))=dQdz(_I3(Nz-1,:,:))
85
86	!?? fluxWest = XA(Q-qXM1)pDdxOp*umask
87	!?? fluxSouth = YA(Q-qYM1)pDdyOp*vmask
88	!?? IF ( tracerBhDiffusion ) fluxUpper=rPVolPMASK(
89	!?? & Cshift(fluxWest,DIM=_I3X,SHIFT=+1)
90	!?? & -fluxWest+Cshift(fluxSouth,DIM=_I3Y,SHIFT=+1)-fluxSouth)
91	IF ( tracerBhDiffusion ) THEN
92	fluxWest = XA(Q-qXM1)pDdxOp*umask
93	fluxSouth = YA(Q-qYM1)pDdyOp*vmask
94	fluxUpper=rPVolPMASK( Cshift(fluxWest,DIM=_I3X,SHIFT=+1)
95	& -fluxWest+Cshift(fluxSouth,DIM=_I3Y,SHIFT=+1)-fluxSouth)
96	ENDIF
97
98	fluxWest = -UMASKXA(
99	& (0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3X,SHIFT=-1))+a2TracerXY)
100	& *dQDx
101	& )
102
103	fluxSouth = -VMASKYA(
104	& (0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3Y,SHIFT=-1))+a2TracerXY)
105	& *dQDy
106	& )
107
108	IF ( tracerBhDiffusion ) THEN
109	fluxWest = fluxWest +
110	& a4TracerXYXA(fluxUpper-Cshift(fluxUpper,DIM=_I3X,SHIFT=-1))*pDdxOp
111	fluxSouth = fluxSouth+
112	& a4TracerXYYA(fluxUpper-Cshift(fluxUpper,DIM=_I3Y,SHIFT=-1))*pDdyOp
113	ENDIF
114
115	! Interpolate d/dx Q and d/dy Q to W points
116	dQdx=0.5d0(dQdx+CSHIFT(dQdx,DIM=_I3X,SHIFT=+1))PMASK ! Interp to T
117	dQdx=0.5d0(dQdx+CSHIFT(dQdx,DIM=_I3Z,SHIFT=-1))WMASK ! Interp to W
118	dQdy=0.5d0(dQdy+CSHIFT(dQdy,DIM=_I3Y,SHIFT=+1))PMASK ! Interp to T
119	dQdy=0.5d0(dQdy+CSHIFT(dQdy,DIM=_I3Z,SHIFT=-1))WMASK ! Interp to W
120
121	!! fluxUpper = -ZA0.5(Kredi+CSHIFT(Kredi,DIM=_I3Z,SHIFT=-1))
122	!! & ( K31dQdx + K32dQdy + K33dQdz )
123
124	fluxUpper = -ZA*(
125	& 0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3Z,SHIFT=-1))
126	& ( K31dQdx + K32*dQdy )
127	& +0.5*(Kredigm+CSHIFT(Kredigm,DIM=_I3Z,SHIFT=-1))
128	& ( K33dQdz )
129	& )
130
131	! Add in vertical diffusion (Av d^2/dz^2)
132	! fluxUpper=-wmaskgronilfluxUpper-a2TracerZwmaskZA(Q-qZM1))*pDdzOP
133	fluxUpper = -WMASK( (GRONIL)fluxUpper-a2TracerZZA(qZM1-Q)pDdzOP )
134	fluxUpper(_I3(1,:,:)) = 0.
135	ENDIF
136	IF ( tracerAdvection ) THEN
137	C Include transport of temperature by flow field.
138	fluxWest = fluxWest +XA(qXM10.5+Q0.5)U*UMASK
139	fluxSouth = fluxSouth+YA(qYM10.5+Q0.5)V*VMASK
140	fluxUpper = fluxUpper+ZA(qZM10.5+Q0.5)W*WMASK
141	ENDIF
142	IF ( tracerDiffusion .OR. tracerAdvection ) THEN
143	qT = qT - (
144	& +Cshift(fluxWest,DIM=_I3X,SHIFT=+1)-fluxWest
145	& +Cshift(fluxSouth,DIM=_I3Y,SHIFT=+1)-fluxSouth
146	& -fluxUpper+Cshift(fluxUpper,DIM=_I3Z,SHIFT=+1)
147	& )*rPVol
148	ENDIF
149	qT(_I3(1,:,:)) = qT(_I3(1,:,:))
150	& + freeSurfFacW(_I3(1,:,:))Q(_I3(1,:,:))/delps(1)
151
152	IF(tracerForcing)THEN
153	C Include external forcing ( e.g. surface heat fluxes, boundary fluxes )
154	C Switch in the appropriate tracer forcing module
155	qT(_I3(1,:,:)) = qT(_I3(1,:,:))
156	& - 1./(30.oneDay)(Q(_I3(1,:,:))-surfProf(:,:))
157	c & - surfProf(_I3(2,:,:))
158	ENDIF
159	qT = qT*PMASK
160
161	C
162	Q = Q + ( (1.5+abEps)qT-(0.5+abEps)qTNM1 )delTasyncfac
163	C
164
165	RETURN
166	END