/[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.5 - (hide annotations) (download)
Fri Feb 3 14:08:28 2012 UTC (13 years, 6 months ago) by dimitri
Branch: MAIN
Changes since 1.4: +6 -2 lines
adding 2-way coupling

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

  ViewVC Help
Powered by ViewVC 1.1.22