/[MITgcm]/MITgcm_contrib/bbl/code/do_oceanic_phys.F
ViewVC logotype

Annotation of /MITgcm_contrib/bbl/code/do_oceanic_phys.F

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


Revision 1.1 - (hide annotations) (download)
Sun Mar 6 02:11:05 2011 UTC (14 years, 4 months ago) by dimitri
Branch: MAIN
changing density comparisons from potential to in situ
code contributed by Madeline Miller

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

  ViewVC Help
Powered by ViewVC 1.1.22