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