/[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.88 - (hide annotations) (download)
Sun Mar 28 14:15:19 2010 UTC (14 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint62f, checkpoint62e
Changes since 1.87: +2 -1 lines
Add one store

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

  ViewVC Help
Powered by ViewVC 1.1.22