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