/[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.83 - (hide annotations) (download)
Tue Dec 15 17:03:28 2009 UTC (14 years, 5 months ago) by jahn
Branch: MAIN
Changes since 1.82: +3 -7 lines
move bi,bj loops into obcs_calc, so obcs_prescribe_read is called only once

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

  ViewVC Help
Powered by ViewVC 1.1.22