/[MITgcm]/MITgcm/model/src/do_oceanic_phys.F
ViewVC logotype

Annotation of /MITgcm/model/src/do_oceanic_phys.F

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

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