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