/[MITgcm]/MITgcm/pkg/ptracers/ptracers_integrate.F
ViewVC logotype

Contents of /MITgcm/pkg/ptracers/ptracers_integrate.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.31 - (show annotations) (download)
Tue Apr 4 14:52:43 2006 UTC (18 years, 1 month ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58f_post, checkpoint58d_post, checkpoint58e_post, checkpoint58g_post, checkpoint58h_post, checkpoint58j_post, checkpoint58i_post
Changes since 1.30: +42 -11 lines
Fix ptracers adjoint
o avoid extensive recomputatations
o fix missing re-init. of gptr (missed by TAF)

1 C $Header: /u/gcmpack/MITgcm/pkg/ptracers/ptracers_integrate.F,v 1.30 2006/03/07 15:28:39 jmc Exp $
2 C $Name: $
3
4 #include "PTRACERS_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: PTRACERS_INTEGRATE
8
9 C !INTERFACE: ==========================================================
10 SUBROUTINE PTRACERS_INTEGRATE(
11 I bi,bj,k,
12 I xA,yA,uTrans,vTrans,rTrans,
13 I rTransKp1,maskUp,
14 X rFlx,KappaRtr,
15 I myIter,myTime,myThid )
16
17 C !DESCRIPTION:
18 C Calculates tendency for passive tracers and integrates forward
19 C in time.
20
21 C !USES: ===============================================================
22 IMPLICIT NONE
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "DYNVARS.h"
27 #include "PTRACERS_SIZE.h"
28 #include "PTRACERS.h"
29 #include "GAD.h"
30 #ifdef ALLOW_AUTODIFF_TAMC
31 # include "tamc.h"
32 # include "tamc_keys.h"
33 #endif
34
35 C !INPUT PARAMETERS: ===================================================
36 C bi,bj :: tile indices
37 C k :: vertical level number
38 C xA :: face area at U points in level k
39 C yA :: face area at V points in level k
40 C uTrans :: zonal transport in level k
41 C vTrans :: meridional transport in level k
42 C rTrans :: vertical volume transport at interface k
43 C rTransKp1 :: vertical volume transport at interface k+1
44 C maskUp :: mask for vertical transport
45 C rFlx :: vertical flux
46 C KappaRtr :: vertical diffusion of passive tracers, interf k
47 C myIter :: time-step number
48 C myTime :: model time
49 C myThid :: thread number
50 INTEGER bi,bj,k
51 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56 _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58 _RL rFlx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
59 _RL KappaRtr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
60 INTEGER myIter
61 _RL myTime
62 INTEGER myThid
63
64 C !OUTPUT PARAMETERS: ==================================================
65 C none
66
67 #ifdef ALLOW_PTRACERS
68
69 C !LOCAL VARIABLES: ====================================================
70 C iTracer :: tracer index
71 C iMin,iMax,jMin,jMax :: loop ranges
72 C kUp,kDown :: toggle indices for even/odd level fluxes
73 C km1 :: =min(1,k-1)
74 C GAD_TR :: passive tracer id (GAD_TR1+iTracer-1)
75 INTEGER iTracer
76 INTEGER iMin,iMax,jMin,jMax
77 INTEGER kUp,kDown,km1
78 INTEGER GAD_TR
79 LOGICAL calcAdvection
80 INTEGER iterNb
81 CEOP
82
83 C Loop ranges for daughter routines
84 iMin = 1-OLx+2
85 iMax = sNx+OLx-1
86 jMin = 1-OLy+2
87 jMax = sNy+OLy-1
88
89 km1 = MAX(1,k-1)
90 kUp = 1+MOD(k+1,2)
91 kDown= 1+MOD(k,2)
92
93 C Loop over tracers
94 DO iTracer=1,PTRACERS_numInUse
95
96 #ifdef ALLOW_AUTODIFF_TAMC
97 act0 = iTracer - 1
98 max0 = PTRACERS_num
99 act1 = bi - myBxLo(myThid)
100 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
101 act2 = bj - myByLo(myThid)
102 max2 = myByHi(myThid) - myByLo(myThid) + 1
103 act3 = myThid - 1
104 max3 = nTx*nTy
105 act4 = ikey_dynamics - 1
106 iptrkey = (act0 + 1)
107 & + act1*max0
108 & + act2*max0*max1
109 & + act3*max0*max1*max2
110 & + act4*max0*max1*max2*max3
111 kkey = (iptrkey-1)*Nr + k
112 #endif /* ALLOW_AUTODIFF_TAMC */
113
114 #ifdef ALLOW_AUTODIFF_TAMC
115 rFlx(1,1,kDown,iTracer) = rFlx(1,1,kDown,iTracer)
116 c
117 CADJ STORE ptracer(:,:,k,bi,bj,iTracer)
118 CADJ & = comlev1_bibj_k_ptracers, key=kkey, byte=isbyte
119 CADJ STORE gPtrnm1(:,:,k,bi,bj,iTracer)
120 CADJ & = comlev1_bibj_k_ptracers, key=kkey, byte=isbyte
121 #endif /* ALLOW_AUTODIFF_TAMC */
122
123 C Calculate active tracer tendencies (gPtr) due to internal processes
124 C (advection, [explicit] diffusion, parameterizations,...)
125 calcAdvection = .NOT.multiDimAdvection
126 & .OR. PTRACERS_advScheme(iTracer).EQ.ENUM_CENTERED_2ND
127 & .OR. PTRACERS_advScheme(iTracer).EQ.ENUM_UPWIND_3RD
128 & .OR. PTRACERS_advScheme(iTracer).EQ.ENUM_CENTERED_4TH
129 GAD_TR = GAD_TR1 + iTracer - 1
130 CALL GAD_CALC_RHS(
131 I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
132 I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
133 I uVel, vVel, wVel,
134 I PTRACERS_diffKh(iTracer),
135 I PTRACERS_diffK4(iTracer),
136 I KappaRtr(1-Olx,1-Oly,iTracer),
137 I gPtrNm1(1-Olx,1-Oly,1,1,1,iTracer),
138 I pTracer(1-Olx,1-Oly,1,1,1,iTracer),
139 I GAD_TR,
140 I PTRACERS_advScheme(iTracer),
141 I PTRACERS_advScheme(iTracer),
142 I calcAdvection, PTRACERS_ImplVertAdv(iTracer),
143 I .FALSE.,
144 U rFlx(1-Olx,1-Oly,1,iTracer),
145 U gPtr(1-Olx,1-Oly,1,1,1,iTracer),
146 I myTime, myIter, myThid )
147
148 C External forcing term(s)
149 IF ( tracForcingOutAB.NE.1 )
150 & CALL PTRACERS_FORCING(
151 I bi,bj,iMin,iMax,jMin,jMax,k,iTracer,
152 U gPtr(1-Olx,1-Oly,1,1,1,iTracer),
153 I surfaceForcingPtr(1-Olx,1-Oly,1,1,iTracer),
154 I myIter,myTime,myThid)
155
156 C If using Adams-Bashforth II, then extrapolate tendencies
157 C gPtr is now the tracer tendency for explicit advection/diffusion
158 IF ( PTRACERS_advScheme(iTracer).EQ.ENUM_CENTERED_2ND
159 & .OR.PTRACERS_advScheme(iTracer).EQ.ENUM_UPWIND_3RD
160 & .OR.PTRACERS_advScheme(iTracer).EQ.ENUM_CENTERED_4TH ) THEN
161 #ifdef ALLOW_MATRIX
162 C If matrix is being computed, block call to S/R ADAMS_BASHFORTH2 to
163 C prevent gPtr from being replaced by the average of gPtr and gPtrNm1.
164 IF (.NOT.useMATRIX) THEN
165 #endif
166 iterNb = myIter
167 IF (staggerTimeStep) iterNb = myIter - 1
168 CALL ADAMS_BASHFORTH2(
169 I bi,bj,K,
170 U gPtr(1-Olx,1-Oly,1,1,1,iTracer),
171 U gPtrNm1(1-Olx,1-Oly,1,1,1,iTracer),
172 I iterNb, myThid )
173 #ifdef ALLOW_MATRIX
174 ENDIF
175 #endif
176 ENDIF
177
178 C External forcing term(s)
179 IF ( tracForcingOutAB.EQ.1 )
180 & CALL PTRACERS_FORCING(
181 I bi,bj,iMin,iMax,jMin,jMax,k,iTracer,
182 U gPtr(1-Olx,1-Oly,1,1,1,iTracer),
183 I surfaceForcingPtr(1-Olx,1-Oly,1,1,iTracer),
184 I myIter,myTime,myThid)
185
186 #ifdef NONLIN_FRSURF
187 C Account for change in level thickness
188 IF (nonlinFreeSurf.GT.0) THEN
189 CALL FREESURF_RESCALE_G(
190 I bi,bj,K,
191 U gPtr(1-Olx,1-Oly,1,1,1,iTracer),
192 I myThid )
193 IF ( PTRACERS_advScheme(iTracer).EQ.ENUM_CENTERED_2ND
194 & .OR.PTRACERS_advScheme(iTracer).EQ.ENUM_UPWIND_3RD
195 & .OR.PTRACERS_advScheme(iTracer).EQ.ENUM_CENTERED_4TH )
196 & CALL FREESURF_RESCALE_G(
197 I bi,bj,K,
198 U gPtrNm1(1-Olx,1-Oly,1,1,1,iTracer),
199 I myThid )
200 ENDIF
201 #endif /* NONLIN_FRSURF */
202
203 C Integrate forward in time, storing in gPtr: G=T+dt*G
204 CALL TIMESTEP_TRACER(
205 I bi,bj,iMin,iMax,jMin,jMax,k,
206 I PTRACERS_advScheme(iTracer),
207 I pTracer(1-Olx,1-Oly,1,1,1,iTracer),
208 I gPtr(1-Olx,1-Oly,1,1,1,iTracer),
209 I myIter,myThid )
210
211 #ifdef ALLOW_OBCS
212 C Apply open boundary conditions
213 IF (useOBCS) THEN
214 CALL OBCS_APPLY_PTRACER(
215 I bi, bj, k, iTracer,
216 U gPtr(1-Olx,1-Oly,k,bi,bj,iTracer),
217 I myThid )
218 END IF
219 #endif /* ALLOW_OBCS */
220 C end of tracer loop
221 ENDDO
222
223 #endif /* ALLOW_PTRACERS */
224
225 RETURN
226 END

  ViewVC Help
Powered by ViewVC 1.1.22