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