/[MITgcm]/MITgcm_contrib/MPMice/beaufort/code/do_oceanic_phys.F
ViewVC logotype

Annotation of /MITgcm_contrib/MPMice/beaufort/code/do_oceanic_phys.F

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

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