/[MITgcm]/MITgcm/model/src/do_oceanic_phys.F
ViewVC logotype

Annotation of /MITgcm/model/src/do_oceanic_phys.F

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


Revision 1.101 - (hide annotations) (download)
Sat Mar 5 17:51:19 2011 UTC (13 years, 3 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint62v, checkpoint62u, checkpoint62t
Changes since 1.100: +16 -3 lines
Some store fixes (probably benign)

1 heimbach 1.101 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.100 2011/03/03 16:53:45 heimbach Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6    
7     #ifdef ALLOW_AUTODIFF_TAMC
8     # ifdef ALLOW_GMREDI
9     # include "GMREDI_OPTIONS.h"
10     # endif
11     # ifdef ALLOW_KPP
12     # include "KPP_OPTIONS.h"
13     # endif
14 jmc 1.29 # ifdef ALLOW_SEAICE
15     # include "SEAICE_OPTIONS.h"
16     # endif
17 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
18    
19     CBOP
20     C !ROUTINE: DO_OCEANIC_PHYS
21     C !INTERFACE:
22     SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
23     C !DESCRIPTION: \bv
24     C *==========================================================*
25 jmc 1.28 C | SUBROUTINE DO_OCEANIC_PHYS
26     C | o Controlling routine for oceanic physics and
27 jmc 1.1 C | parameterization
28     C *==========================================================*
29     C | o originally, part of S/R thermodynamics
30     C *==========================================================*
31     C \ev
32    
33     C !USES:
34     IMPLICIT NONE
35     C == Global variables ===
36     #include "SIZE.h"
37     #include "EEPARAMS.h"
38     #include "PARAMS.h"
39 jmc 1.69 #include "GRID.h"
40 jmc 1.1 #include "DYNVARS.h"
41 jmc 1.20 #ifdef ALLOW_TIMEAVE
42     #include "TIMEAVE_STATV.h"
43     #endif
44 mlosch 1.22 #if defined (ALLOW_BALANCE_FLUXES) && !(defined ALLOW_AUTODIFF_TAMC)
45     #include "FFIELDS.h"
46     #endif
47 jmc 1.1
48     #ifdef ALLOW_AUTODIFF_TAMC
49 heimbach 1.80 # include "AUTODIFF_MYFIELDS.h"
50 jmc 1.1 # include "tamc.h"
51     # include "tamc_keys.h"
52     # include "FFIELDS.h"
53 heimbach 1.54 # include "SURFACE.h"
54 jmc 1.1 # include "EOS.h"
55     # ifdef ALLOW_KPP
56     # include "KPP.h"
57     # endif
58 gforget 1.91 # ifdef ALLOW_GGL90
59     # include "GGL90.h"
60     # endif
61 jmc 1.1 # ifdef ALLOW_GMREDI
62     # include "GMREDI.h"
63     # endif
64     # ifdef ALLOW_EBM
65     # include "EBM.h"
66     # endif
67 jmc 1.29 # ifdef ALLOW_EXF
68     # include "ctrl.h"
69 jmc 1.40 # include "EXF_FIELDS.h"
70 jmc 1.29 # ifdef ALLOW_BULKFORMULAE
71 jmc 1.40 # include "EXF_CONSTANTS.h"
72 jmc 1.29 # endif
73     # endif
74     # ifdef ALLOW_SEAICE
75     # include "SEAICE.h"
76     # endif
77 heimbach 1.75 # ifdef ALLOW_SALT_PLUME
78     # include "SALT_PLUME.h"
79     # endif
80 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
81    
82     C !INPUT/OUTPUT PARAMETERS:
83     C == Routine arguments ==
84 jmc 1.18 C myTime :: Current time in simulation
85     C myIter :: Current iteration number in simulation
86     C myThid :: Thread number for this instance of the routine.
87 jmc 1.1 _RL myTime
88     INTEGER myIter
89     INTEGER myThid
90    
91     C !LOCAL VARIABLES:
92     C == Local variables
93 jmc 1.47 C rhoK, rhoKm1 :: Density at current level, and level above
94 jmc 1.18 C iMin, iMax :: Ranges and sub-block indices on which calculations
95 jmc 1.1 C jMin, jMax are applied.
96 jmc 1.18 C bi, bj :: tile indices
97     C i,j,k :: loop indices
98 jmc 1.47 _RL rhoKp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 jmc 1.1 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103     INTEGER iMin, iMax
104     INTEGER jMin, jMax
105     INTEGER bi, bj
106 jmc 1.18 INTEGER i, j, k
107 jmc 1.17 INTEGER doDiagsRho
108     #ifdef ALLOW_DIAGNOSTICS
109     LOGICAL DIAGNOSTICS_IS_ON
110     EXTERNAL DIAGNOSTICS_IS_ON
111     #endif /* ALLOW_DIAGNOSTICS */
112 jmc 1.1
113     CEOP
114    
115 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
116     C-- dummy statement to end declaration part
117     itdkey = 1
118     #endif /* ALLOW_AUTODIFF_TAMC */
119    
120 jmc 1.1 #ifdef ALLOW_DEBUG
121 jmc 1.96 IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
122 jmc 1.1 #endif
123 jmc 1.36
124 jmc 1.17 doDiagsRho = 0
125     #ifdef ALLOW_DIAGNOSTICS
126     IF ( useDiagnostics .AND. fluidIsWater ) THEN
127 jmc 1.74 IF ( DIAGNOSTICS_IS_ON('WRHOMASS',myThid) )
128 jmc 1.71 & doDiagsRho = doDiagsRho + 1
129     IF ( DIAGNOSTICS_IS_ON('DRHODR ',myThid) )
130     & doDiagsRho = doDiagsRho + 2
131     IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
132     & doDiagsRho = doDiagsRho + 4
133 jmc 1.17 ENDIF
134     #endif /* ALLOW_DIAGNOSTICS */
135    
136 jmc 1.82 #ifdef ALLOW_OBCS
137     IF (useOBCS) THEN
138     C-- Calculate future values on open boundaries
139     C-- moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
140 heimbach 1.100 # ifdef ALLOW_AUTODIFF_TAMC
141     CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
142     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
143     # endif
144     # ifdef ALLOW_DEBUG
145 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
146 heimbach 1.100 # endif
147 jahn 1.83 CALL OBCS_CALC( myTime+deltaTclock, myIter+1,
148     I uVel, vVel, wVel, theta, salt, myThid )
149 jmc 1.82 ENDIF
150     #endif /* ALLOW_OBCS */
151 jmc 1.69
152 jmc 1.98 #ifdef ALLOW_ADDFLUID
153     IF ( fluidIsWater ) THEN
154     DO bj=myByLo(myThid),myByHi(myThid)
155     DO bi=myBxLo(myThid),myBxHi(myThid)
156     DO k=1,Nr
157     DO j=1-OLy,sNy+OLy
158     DO i=1-OLx,sNx+OLx
159     addMass(i,j,k,bi,bj) = 0. _d 0
160     ENDDO
161     ENDDO
162     ENDDO
163     ENDDO
164     ENDDO
165     ENDIF
166     #endif /* ALLOW_ADDFLUID */
167    
168 gforget 1.87 #ifdef ALLOW_AUTODIFF_TAMC
169     # ifdef ALLOW_SALT_PLUME
170     DO bj=myByLo(myThid),myByHi(myThid)
171     DO bi=myBxLo(myThid),myBxHi(myThid)
172     DO j=1-OLy,sNy+OLy
173     DO i=1-OLx,sNx+OLx
174     saltPlumeDepth(i,j,bi,bj) = 0. _d 0
175     saltPlumeFlux(i,j,bi,bj) = 0. _d 0
176     ENDDO
177     ENDDO
178     ENDDO
179     ENDDO
180     # endif
181     #endif /* ALLOW_AUTODIFF_TAMC */
182    
183 jmc 1.29 #ifdef ALLOW_SEAICE
184     IF ( useSEAICE ) THEN
185 heimbach 1.62 # ifdef ALLOW_AUTODIFF_TAMC
186 heimbach 1.65 cph-adj-test(
187 heimbach 1.81 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
188     CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
189 heimbach 1.88 CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte
190 heimbach 1.81 CADJ STORE empmr,qsw,theta = comlev1, key = ikey_dynamics,
191 heimbach 1.77 CADJ & kind = isbyte
192 heimbach 1.65 cph-adj-test)
193 heimbach 1.77 CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
194     CADJ & kind = isbyte
195     CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
196     CADJ & kind = isbyte
197 heimbach 1.34 cph# ifdef EXF_READ_EVAP
198 heimbach 1.77 CADJ STORE evap = comlev1, key = ikey_dynamics,
199     CADJ & kind = isbyte
200 heimbach 1.34 cph# endif
201 heimbach 1.77 CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
202     CADJ & kind = isbyte
203 heimbach 1.95 # ifdef SEAICE_CGRID
204 heimbach 1.93 CADJ STORE stressdivergencex = comlev1, key = ikey_dynamics,
205     CADJ & kind = isbyte
206     CADJ STORE stressdivergencey = comlev1, key = ikey_dynamics,
207     CADJ & kind = isbyte
208     # endif
209 heimbach 1.62 # ifdef SEAICE_ALLOW_DYNAMICS
210 heimbach 1.77 CADJ STORE uice = comlev1, key = ikey_dynamics,
211     CADJ & kind = isbyte
212     CADJ STORE vice = comlev1, key = ikey_dynamics,
213     CADJ & kind = isbyte
214 heimbach 1.62 # ifdef SEAICE_ALLOW_EVP
215 heimbach 1.77 CADJ STORE seaice_sigma1 = comlev1, key = ikey_dynamics,
216     CADJ & kind = isbyte
217     CADJ STORE seaice_sigma2 = comlev1, key = ikey_dynamics,
218     CADJ & kind = isbyte
219     CADJ STORE seaice_sigma12 = comlev1, key = ikey_dynamics,
220     CADJ & kind = isbyte
221 heimbach 1.62 # endif
222     # endif
223     # ifdef SEAICE_SALINITY
224 heimbach 1.77 CADJ STORE salt = comlev1, key = ikey_dynamics,
225     CADJ & kind = isbyte
226 heimbach 1.62 # endif
227     # ifdef ATMOSPHERIC_LOADING
228 heimbach 1.77 CADJ STORE pload = comlev1, key = ikey_dynamics,
229     CADJ & kind = isbyte
230     CADJ STORE siceload = comlev1, key = ikey_dynamics,
231     CADJ & kind = isbyte
232 heimbach 1.62 # endif
233     # ifdef NONLIN_FRSURF
234 heimbach 1.77 CADJ STORE recip_hfacc = comlev1, key = ikey_dynamics,
235     CADJ & kind = isbyte
236 heimbach 1.62 # endif
237 heimbach 1.78 # ifdef ANNUAL_BALANCE
238     CADJ STORE balance_itcount = comlev1, key = ikey_dynamics,
239     CADJ & kind = isbyte
240     # endif /* ANNUAL_BALANCE */
241 heimbach 1.55 # endif
242 heimbach 1.62 # ifdef ALLOW_DEBUG
243 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
244 heimbach 1.62 # endif
245 jmc 1.29 CALL TIMER_START('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
246     CALL SEAICE_MODEL( myTime, myIter, myThid )
247     CALL TIMER_STOP ('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
248 heimbach 1.62 # ifdef ALLOW_COST
249 heimbach 1.57 CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
250 heimbach 1.62 # endif
251 heimbach 1.35 ENDIF
252 jmc 1.29 #endif /* ALLOW_SEAICE */
253    
254 heimbach 1.64 #ifdef ALLOW_AUTODIFF_TAMC
255 heimbach 1.77 CADJ STORE sst, sss = comlev1, key = ikey_dynamics,
256     CADJ & kind = isbyte
257     CADJ STORE qsw = comlev1, key = ikey_dynamics,
258     CADJ & kind = isbyte
259 heimbach 1.64 # ifdef ALLOW_SEAICE
260 heimbach 1.77 CADJ STORE area = comlev1, key = ikey_dynamics,
261     CADJ & kind = isbyte
262 heimbach 1.64 # endif
263     #endif
264    
265 jscott 1.30 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
266 jmc 1.14 IF ( useThSIce .AND. fluidIsWater ) THEN
267 heimbach 1.101 # ifdef ALLOW_AUTODIFF_TAMC
268     cph(
269     # ifdef NONLIN_FRSURF
270     CADJ STORE uice,vice = comlev1, key = ikey_dynamics,
271     CADJ & kind = isbyte
272     CADJ STORE salt,theta = comlev1, key = ikey_dynamics,
273     CADJ & kind = isbyte
274     CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics,
275     CADJ & kind = isbyte
276     CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics,
277     CADJ & kind = isbyte
278     # endif
279     # endif
280     # ifdef ALLOW_DEBUG
281 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
282 heimbach 1.101 # endif
283 jmc 1.5 C-- Step forward Therm.Sea-Ice variables
284     C and modify forcing terms including effects from ice
285     CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
286     CALL THSICE_MAIN( myTime, myIter, myThid )
287     CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
288     ENDIF
289     #endif /* ALLOW_THSICE */
290    
291 mlosch 1.21 #ifdef ALLOW_SHELFICE
292 heimbach 1.92 # ifdef ALLOW_AUTODIFF_TAMC
293     CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
294     CADJ & kind = isbyte
295     # endif
296 mlosch 1.21 IF ( useShelfIce .AND. fluidIsWater ) THEN
297     #ifdef ALLOW_DEBUG
298 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
299 mlosch 1.21 #endif
300 jmc 1.47 C compute temperature and (virtual) salt flux at the
301 mlosch 1.21 C shelf-ice ocean interface
302     CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
303     & myThid)
304     CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
305     CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
306     & myThid)
307     ENDIF
308     #endif /* ALLOW_SHELFICE */
309    
310 dimitri 1.85 #ifdef ALLOW_ICEFRONT
311     IF ( useICEFRONT .AND. fluidIsWater ) THEN
312     #ifdef ALLOW_DEBUG
313 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
314 dimitri 1.85 #endif
315     C compute temperature and (virtual) salt flux at the
316     C ice-front ocean interface
317     CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
318     & myThid)
319     CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
320     CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
321     & myThid)
322     ENDIF
323     #endif /* ALLOW_ICEFRONT */
324    
325 jmc 1.5 C-- Freeze water at the surface
326     #ifdef ALLOW_AUTODIFF_TAMC
327 heimbach 1.77 CADJ STORE theta = comlev1, key = ikey_dynamics,
328     CADJ & kind = isbyte
329 jmc 1.5 #endif
330 jmc 1.70 IF ( allowFreezing ) THEN
331 jmc 1.5 CALL FREEZE_SURFACE( myTime, myIter, myThid )
332     ENDIF
333    
334 jmc 1.28 #ifdef ALLOW_OCN_COMPON_INTERF
335 jmc 1.5 C-- Apply imported data (from coupled interface) to forcing fields
336 jmc 1.28 C jmc: do not know precisely where to put this call (bf or af thSIce ?)
337 jmc 1.36 IF ( useCoupler ) THEN
338 jmc 1.5 CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
339 jmc 1.36 ENDIF
340 jmc 1.28 #endif /* ALLOW_OCN_COMPON_INTERF */
341 jmc 1.5
342 jmc 1.25 #ifdef ALLOW_BALANCE_FLUXES
343 jmc 1.36 C balance fluxes
344     IF ( balanceEmPmR )
345 jmc 1.84 & CALL REMOVE_MEAN_RS( 1, EmPmR, maskInC, maskInC, rA, drF,
346 jmc 1.25 & 'EmPmR', myTime, myThid )
347 jmc 1.36 IF ( balanceQnet )
348 jmc 1.84 & CALL REMOVE_MEAN_RS( 1, Qnet, maskInC, maskInC, rA, drF,
349 jmc 1.25 & 'Qnet ', myTime, myThid )
350     #endif /* ALLOW_BALANCE_FLUXES */
351    
352 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
353     C-- HPF directive to help TAMC
354     CHPF$ INDEPENDENT
355     #endif /* ALLOW_AUTODIFF_TAMC */
356     DO bj=myByLo(myThid),myByHi(myThid)
357     #ifdef ALLOW_AUTODIFF_TAMC
358 heimbach 1.15 C-- HPF directive to help TAMC
359     CHPF$ INDEPENDENT
360 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
361     DO bi=myBxLo(myThid),myBxHi(myThid)
362    
363     #ifdef ALLOW_AUTODIFF_TAMC
364     act1 = bi - myBxLo(myThid)
365     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
366     act2 = bj - myByLo(myThid)
367     max2 = myByHi(myThid) - myByLo(myThid) + 1
368     act3 = myThid - 1
369     max3 = nTx*nTy
370     act4 = ikey_dynamics - 1
371     itdkey = (act1 + 1) + act2*max1
372     & + act3*max1*max2
373     & + act4*max1*max2*max3
374 jmc 1.74 #else /* ALLOW_AUTODIFF_TAMC */
375 jmc 1.47 C if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
376 jmc 1.36 C and all vertical mixing schemes, but keep OBCS_CALC
377     IF ( fluidIsWater ) THEN
378 jmc 1.74 #endif /* ALLOW_AUTODIFF_TAMC */
379 jmc 1.1
380     C-- Set up work arrays with valid (i.e. not NaN) values
381     C These inital values do not alter the numerical results. They
382     C just ensure that all memory references are to valid floating
383     C point numbers. This prevents spurious hardware signals due to
384     C uninitialised but inert locations.
385    
386 jmc 1.69 #ifdef ALLOW_AUTODIFF_TAMC
387 jmc 1.1 DO j=1-OLy,sNy+OLy
388     DO i=1-OLx,sNx+OLx
389 jmc 1.69 rhoKm1 (i,j) = 0. _d 0
390 jmc 1.47 rhoKp1 (i,j) = 0. _d 0
391 jmc 1.1 ENDDO
392     ENDDO
393 jmc 1.69 #endif /* ALLOW_AUTODIFF_TAMC */
394 jmc 1.1
395     DO k=1,Nr
396     DO j=1-OLy,sNy+OLy
397     DO i=1-OLx,sNx+OLx
398 jmc 1.69 C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
399 jmc 1.1 sigmaX(i,j,k) = 0. _d 0
400     sigmaY(i,j,k) = 0. _d 0
401     sigmaR(i,j,k) = 0. _d 0
402     #ifdef ALLOW_AUTODIFF_TAMC
403     cph all the following init. are necessary for TAF
404     cph although some of these are re-initialised later.
405 jmc 1.71 c rhoInSitu(i,j,k,bi,bj) = 0.
406 jmc 1.1 IVDConvCount(i,j,k,bi,bj) = 0.
407     # ifdef ALLOW_GMREDI
408     Kwx(i,j,k,bi,bj) = 0. _d 0
409     Kwy(i,j,k,bi,bj) = 0. _d 0
410     Kwz(i,j,k,bi,bj) = 0. _d 0
411     # ifdef GM_NON_UNITY_DIAGONAL
412     Kux(i,j,k,bi,bj) = 0. _d 0
413     Kvy(i,j,k,bi,bj) = 0. _d 0
414     # endif
415     # ifdef GM_EXTRA_DIAGONAL
416     Kuz(i,j,k,bi,bj) = 0. _d 0
417     Kvz(i,j,k,bi,bj) = 0. _d 0
418     # endif
419     # ifdef GM_BOLUS_ADVEC
420     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
421     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
422     # endif
423     # ifdef GM_VISBECK_VARIABLE_K
424     VisbeckK(i,j,bi,bj) = 0. _d 0
425     # endif
426     # endif /* ALLOW_GMREDI */
427 heimbach 1.42 # ifdef ALLOW_KPP
428     KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
429     KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
430     # endif /* ALLOW_KPP */
431 gforget 1.91 # ifdef ALLOW_GGL90
432     GGL90viscArU(i,j,k,bi,bj) = 0. _d 0
433     GGL90viscArV(i,j,k,bi,bj) = 0. _d 0
434     GGL90diffKr(i,j,k,bi,bj) = 0. _d 0
435     # endif /* ALLOW_GGL90 */
436 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
437     ENDDO
438     ENDDO
439     ENDDO
440    
441     iMin = 1-OLx
442     iMax = sNx+OLx
443     jMin = 1-OLy
444     jMax = sNy+OLy
445    
446     #ifdef ALLOW_AUTODIFF_TAMC
447 jmc 1.96 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
448 heimbach 1.77 CADJ & kind = isbyte
449 jmc 1.96 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
450 heimbach 1.77 CADJ & kind = isbyte
451 heimbach 1.12 CADJ STORE totphihyd(:,:,:,bi,bj)
452 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
453 heimbach 1.77 CADJ & kind = isbyte
454 heimbach 1.10 # ifdef ALLOW_KPP
455 jmc 1.96 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
456 heimbach 1.77 CADJ & kind = isbyte
457 jmc 1.96 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
458 heimbach 1.77 CADJ & kind = isbyte
459 heimbach 1.10 # endif
460 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
461    
462     #ifdef ALLOW_DEBUG
463 jmc 1.96 IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
464 jmc 1.1 #endif
465    
466     C-- Start of diagnostic loop
467     DO k=Nr,1,-1
468    
469     #ifdef ALLOW_AUTODIFF_TAMC
470     C? Patrick, is this formula correct now that we change the loop range?
471     C? Do we still need this?
472 jmc 1.69 cph kkey formula corrected.
473 jmc 1.47 cph Needed for rhoK, rhoKm1, in the case useGMREDI.
474 jmc 1.71 kkey = (itdkey-1)*Nr + k
475 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
476    
477 jmc 1.71 C-- Always compute density (stored in common block) here; even when it is not
478     C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
479     #ifdef ALLOW_DEBUG
480 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D',myThid)
481 jmc 1.71 #endif
482     #ifdef ALLOW_AUTODIFF_TAMC
483 jmc 1.74 IF ( fluidIsWater ) THEN
484 jmc 1.96 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
485 heimbach 1.77 CADJ & kind = isbyte
486 jmc 1.96 CADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
487 heimbach 1.77 CADJ & kind = isbyte
488 jmc 1.71 #endif /* ALLOW_AUTODIFF_TAMC */
489 jmc 1.69 #ifdef ALLOW_DOWN_SLOPE
490     IF ( useDOWN_SLOPE ) THEN
491     CALL DWNSLP_CALC_RHO(
492     I theta, salt,
493 jmc 1.71 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
494 jmc 1.69 I k, bi, bj, myTime, myIter, myThid )
495 jmc 1.71 ELSE
496     #endif /* ALLOW_DOWN_SLOPE */
497     CALL FIND_RHO_2D(
498     I iMin, iMax, jMin, jMax, k,
499     I theta(1-OLx,1-OLy,k,bi,bj),
500     I salt (1-OLx,1-OLy,k,bi,bj),
501     O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
502     I k, bi, bj, myThid )
503     #ifdef ALLOW_DOWN_SLOPE
504 jmc 1.69 ENDIF
505     #endif /* ALLOW_DOWN_SLOPE */
506 jmc 1.74 #ifdef ALLOW_AUTODIFF_TAMC
507     ELSE
508     C- fluid is not water:
509     DO j=1-OLy,sNy+OLy
510     DO i=1-OLx,sNx+OLx
511     rhoInSitu(i,j,k,bi,bj) = 0.
512     ENDDO
513     ENDDO
514     ENDIF
515     #endif /* ALLOW_AUTODIFF_TAMC */
516 jmc 1.69
517 jmc 1.1 C-- Calculate gradients of potential density for isoneutral
518     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
519 jmc 1.17 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
520 dimitri 1.61 & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
521 jmc 1.1 IF (k.GT.1) THEN
522     #ifdef ALLOW_AUTODIFF_TAMC
523 jmc 1.96 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
524 heimbach 1.77 CADJ & kind = isbyte
525 jmc 1.96 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
526 heimbach 1.77 CADJ & kind = isbyte
527 jmc 1.96 CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey,
528 heimbach 1.77 CADJ & kind = isbyte
529 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
530 jmc 1.68 CALL FIND_RHO_2D(
531     I iMin, iMax, jMin, jMax, k,
532     I theta(1-OLx,1-OLy,k-1,bi,bj),
533     I salt (1-OLx,1-OLy,k-1,bi,bj),
534     O rhoKm1,
535     I k-1, bi, bj, myThid )
536 jmc 1.1 ENDIF
537     #ifdef ALLOW_DEBUG
538 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
539 jmc 1.1 #endif
540 heimbach 1.31 cph Avoid variable aliasing for adjoint !!!
541     DO j=jMin,jMax
542     DO i=iMin,iMax
543 jmc 1.71 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
544 heimbach 1.31 ENDDO
545     ENDDO
546 jmc 1.1 CALL GRAD_SIGMA(
547     I bi, bj, iMin, iMax, jMin, jMax, k,
548 jmc 1.71 I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
549 jmc 1.1 O sigmaX, sigmaY, sigmaR,
550     I myThid )
551 gforget 1.66 #ifdef ALLOW_AUTODIFF_TAMC
552 jmc 1.69 #ifdef GMREDI_WITH_STABLE_ADJOINT
553 gforget 1.66 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
554     cgf -> cuts adjoint dependency from slope to state
555 jmc 1.69 CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
556     CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
557     CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
558 gforget 1.66 #endif
559     #endif /* ALLOW_AUTODIFF_TAMC */
560 jmc 1.1 ENDIF
561    
562     C-- Implicit Vertical Diffusion for Convection
563     c ==> should use sigmaR !!!
564     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
565     #ifdef ALLOW_DEBUG
566 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
567 jmc 1.1 #endif
568     CALL CALC_IVDC(
569     I bi, bj, iMin, iMax, jMin, jMax, k,
570 jmc 1.71 I rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
571 jmc 1.1 I myTime, myIter, myThid)
572     ENDIF
573    
574 jmc 1.17 #ifdef ALLOW_DIAGNOSTICS
575 jmc 1.71 IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
576     CALL DIAGS_RHO_L( k, bi, bj,
577     I rhoInSitu(1-OLx,1-OLy,k,bi,bj),
578 jmc 1.74 I rhoKm1, wVel,
579 jmc 1.71 I myTime, myIter, myThid )
580 jmc 1.17 ENDIF
581     #endif
582    
583 jmc 1.1 C-- end of diagnostic k loop (Nr:1)
584     ENDDO
585    
586 heimbach 1.57 #ifdef ALLOW_AUTODIFF_TAMC
587 jmc 1.69 CADJ STORE IVDConvCount(:,:,:,bi,bj)
588 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
589 heimbach 1.77 CADJ & kind = isbyte
590 heimbach 1.57 #endif
591    
592 jmc 1.47 C-- Diagnose Mixed Layer Depth:
593 jmc 1.71 IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
594     CALL CALC_OCE_MXLAYER(
595     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
596     I bi, bj, myTime, myIter, myThid )
597 jmc 1.47 ENDIF
598 heimbach 1.53
599 dimitri 1.52 #ifdef ALLOW_SALT_PLUME
600 dimitri 1.61 IF ( useSALT_PLUME ) THEN
601 jmc 1.71 CALL SALT_PLUME_CALC_DEPTH(
602     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
603     I bi, bj, myTime, myIter, myThid )
604 dimitri 1.60 ENDIF
605 dimitri 1.61 #endif /* ALLOW_SALT_PLUME */
606    
607 jmc 1.8 #ifdef ALLOW_DIAGNOSTICS
608 jmc 1.71 IF ( MOD(doDiagsRho,4).GE.2 ) THEN
609 jmc 1.16 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
610     & 2, bi, bj, myThid)
611 jmc 1.8 ENDIF
612 dimitri 1.61 #endif /* ALLOW_DIAGNOSTICS */
613 jmc 1.8
614 jmc 1.1 C-- Determines forcing terms based on external fields
615     C relaxation terms, etc.
616     #ifdef ALLOW_DEBUG
617 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
618 jmc 1.1 #endif
619 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
620     CADJ STORE EmPmR(:,:,bi,bj)
621 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
622 heimbach 1.77 CADJ & kind = isbyte
623 heimbach 1.26 # ifdef EXACT_CONSERV
624 heimbach 1.23 CADJ STORE PmEpR(:,:,bi,bj)
625 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
626 heimbach 1.77 CADJ & kind = isbyte
627 heimbach 1.26 # endif
628 heimbach 1.27 # ifdef NONLIN_FRSURF
629     CADJ STORE hFac_surfC(:,:,bi,bj)
630 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
631 heimbach 1.77 CADJ & kind = isbyte
632 heimbach 1.27 CADJ STORE recip_hFacC(:,:,:,bi,bj)
633 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
634 heimbach 1.77 CADJ & kind = isbyte
635 heimbach 1.97 # if (defined (ALLOW_PTRACERS))
636     CADJ STORE surfaceForcingS(:,:,bi,bj) = comlev1_bibj, key=itdkey,
637     CADJ & kind = isbyte
638     CADJ STORE surfaceForcingT(:,:,bi,bj) = comlev1_bibj, key=itdkey,
639     CADJ & kind = isbyte
640     # endif /* ALLOW_PTRACERS */
641     # endif /* NONLIN_FRSURF */
642 heimbach 1.23 #endif
643 jmc 1.36 CALL EXTERNAL_FORCING_SURF(
644 jmc 1.1 I bi, bj, iMin, iMax, jMin, jMax,
645     I myTime, myIter, myThid )
646 heimbach 1.27 #ifdef ALLOW_AUTODIFF_TAMC
647     # ifdef EXACT_CONSERV
648     cph-test
649     cphCADJ STORE PmEpR(:,:,bi,bj)
650 jmc 1.96 cphCADJ & = comlev1_bibj, key=itdkey,
651 heimbach 1.77 cphCADJ & kind = isbyte
652 heimbach 1.27 # endif
653     #endif
654 jmc 1.1
655     #ifdef ALLOW_AUTODIFF_TAMC
656     cph needed for KPP
657 jmc 1.4 CADJ STORE surfaceForcingU(:,:,bi,bj)
658 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
659 heimbach 1.77 CADJ & kind = isbyte
660 jmc 1.4 CADJ STORE surfaceForcingV(:,:,bi,bj)
661 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
662 heimbach 1.77 CADJ & kind = isbyte
663 jmc 1.4 CADJ STORE surfaceForcingS(:,:,bi,bj)
664 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
665 heimbach 1.77 CADJ & kind = isbyte
666 jmc 1.4 CADJ STORE surfaceForcingT(:,:,bi,bj)
667 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
668 heimbach 1.77 CADJ & kind = isbyte
669 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
670 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
671 heimbach 1.77 CADJ & kind = isbyte
672 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
673    
674     #ifdef ALLOW_KPP
675     C-- Compute KPP mixing coefficients
676     IF (useKPP) THEN
677     #ifdef ALLOW_DEBUG
678 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
679 jmc 1.1 #endif
680 dfer 1.76 CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
681 jmc 1.1 CALL KPP_CALC(
682 jmc 1.44 I bi, bj, myTime, myIter, myThid )
683 dfer 1.76 CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
684 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
685     ELSE
686     CALL KPP_CALC_DUMMY(
687 jmc 1.44 I bi, bj, myTime, myIter, myThid )
688 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
689     ENDIF
690    
691     #endif /* ALLOW_KPP */
692    
693 mlosch 1.6 #ifdef ALLOW_PP81
694     C-- Compute PP81 mixing coefficients
695     IF (usePP81) THEN
696     #ifdef ALLOW_DEBUG
697 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
698 mlosch 1.6 #endif
699     CALL PP81_CALC(
700     I bi, bj, myTime, myThid )
701     ENDIF
702     #endif /* ALLOW_PP81 */
703    
704     #ifdef ALLOW_MY82
705     C-- Compute MY82 mixing coefficients
706     IF (useMY82) THEN
707     #ifdef ALLOW_DEBUG
708 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
709 mlosch 1.6 #endif
710     CALL MY82_CALC(
711     I bi, bj, myTime, myThid )
712     ENDIF
713     #endif /* ALLOW_MY82 */
714    
715 mlosch 1.9 #ifdef ALLOW_GGL90
716 gforget 1.91 #ifdef ALLOW_AUTODIFF_TAMC
717     CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
718     CADJ & kind = isbyte
719     #endif /* ALLOW_AUTODIFF_TAMC */
720 mlosch 1.9 C-- Compute GGL90 mixing coefficients
721     IF (useGGL90) THEN
722     #ifdef ALLOW_DEBUG
723 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
724 mlosch 1.9 #endif
725 dfer 1.76 CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
726 mlosch 1.9 CALL GGL90_CALC(
727 jmc 1.90 I bi, bj, myTime, myIter, myThid )
728 dfer 1.76 CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
729 mlosch 1.9 ENDIF
730     #endif /* ALLOW_GGL90 */
731    
732 jmc 1.20 #ifdef ALLOW_TIMEAVE
733 jmc 1.36 IF ( taveFreq.GT. 0. _d 0 ) THEN
734 jmc 1.20 CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
735     ENDIF
736     IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
737     CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
738     I Nr, deltaTclock, bi, bj, myThid)
739     ENDIF
740     #endif /* ALLOW_TIMEAVE */
741    
742 jmc 1.69 #ifdef ALLOW_GMREDI
743 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
744     # ifndef GM_EXCLUDE_CLIPPING
745     cph storing here is needed only for one GMREDI_OPTIONS:
746     cph define GM_BOLUS_ADVEC
747     cph keep it although TAF says you dont need to.
748 jmc 1.86 cph but I have avoided the #ifdef for now, in case more things change
749 jmc 1.96 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey,
750 heimbach 1.77 CADJ & kind = isbyte
751 jmc 1.96 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey,
752 heimbach 1.77 CADJ & kind = isbyte
753 jmc 1.96 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey,
754 heimbach 1.77 CADJ & kind = isbyte
755 jmc 1.47 # endif
756     #endif /* ALLOW_AUTODIFF_TAMC */
757    
758     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
759     IF (useGMRedi) THEN
760     #ifdef ALLOW_DEBUG
761 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
762 jmc 1.47 #endif
763     CALL GMREDI_CALC_TENSOR(
764 jmc 1.51 I iMin, iMax, jMin, jMax,
765 jmc 1.47 I sigmaX, sigmaY, sigmaR,
766 jmc 1.51 I bi, bj, myTime, myIter, myThid )
767 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
768     ELSE
769     CALL GMREDI_CALC_TENSOR_DUMMY(
770 jmc 1.51 I iMin, iMax, jMin, jMax,
771 jmc 1.47 I sigmaX, sigmaY, sigmaR,
772 jmc 1.51 I bi, bj, myTime, myIter, myThid )
773 jmc 1.47 #endif /* ALLOW_AUTODIFF_TAMC */
774     ENDIF
775 jmc 1.69 #endif /* ALLOW_GMREDI */
776    
777     #ifdef ALLOW_DOWN_SLOPE
778     IF ( useDOWN_SLOPE ) THEN
779     C-- Calculate Downsloping Flow for Down_Slope parameterization
780     IF ( usingPCoords ) THEN
781     CALL DWNSLP_CALC_FLOW(
782 jmc 1.71 I bi, bj, kSurfC, rhoInSitu,
783 jmc 1.69 I myTime, myIter, myThid )
784     ELSE
785     CALL DWNSLP_CALC_FLOW(
786 jmc 1.71 I bi, bj, kLowC, rhoInSitu,
787 jmc 1.69 I myTime, myIter, myThid )
788     ENDIF
789     ENDIF
790     #endif /* ALLOW_DOWN_SLOPE */
791 jmc 1.47
792 jmc 1.98 #ifndef ALLOW_AUTODIFF_TAMC
793     C--- if fluid Is Water: end
794     ENDIF
795     #endif
796    
797 dimitri 1.94 #ifdef ALLOW_MYPACKAGE
798     IF ( useMYPACKAGE ) THEN
799     CALL MYPACKAGE_CALC_RHS(
800     I bi, bj, myTime, myIter, myThid )
801     ENDIF
802     #endif /* ALLOW_MYPACKAGE */
803    
804 jmc 1.1 C-- end bi,bj loops.
805     ENDDO
806     ENDDO
807    
808 jmc 1.99 #ifdef ALLOW_GMREDI
809     IF ( useGMRedi ) THEN
810     CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
811     ENDIF
812     #endif /* ALLOW_GMREDI */
813    
814 jmc 1.98 #ifdef ALLOW_KPP
815 jmc 1.45 IF (useKPP) THEN
816     CALL KPP_DO_EXCH( myThid )
817     ENDIF
818 jmc 1.98 #endif /* ALLOW_KPP */
819 jmc 1.45
820 jmc 1.18 #ifdef ALLOW_DIAGNOSTICS
821     IF ( fluidIsWater .AND. useDiagnostics ) THEN
822 jmc 1.74 CALL DIAGS_RHO_G(
823 jmc 1.71 I rhoInSitu, uVel, vVel,
824     I myTime, myIter, myThid )
825 jmc 1.18 CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
826     ENDIF
827 jmc 1.19 IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
828 jmc 1.71 CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
829     & 0, Nr, 0, 1, 1, myThid )
830 jmc 1.19 ENDIF
831 jmc 1.18 #endif
832    
833 jmc 1.1 #ifdef ALLOW_DEBUG
834 jmc 1.96 IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
835 jmc 1.1 #endif
836    
837     RETURN
838     END

  ViewVC Help
Powered by ViewVC 1.1.22