/[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.148 - (hide annotations) (download)
Fri Aug 12 14:15:58 2016 UTC (7 years, 9 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint66c, checkpoint66b, checkpoint66a, checkpoint65z
Changes since 1.147: +9 -5 lines
Adjust several STOREs

1 heimbach 1.148 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.147 2015/11/12 20:05:21 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6 jmc 1.134 #ifdef ALLOW_AUTODIFF
7     # include "AUTODIFF_OPTIONS.h"
8     #endif
9     #ifdef ALLOW_CTRL
10     # include "CTRL_OPTIONS.h"
11     #endif
12 heimbach 1.136 #ifdef ALLOW_SALT_PLUME
13     # include "SALT_PLUME_OPTIONS.h"
14     #endif
15 atn 1.142 #ifdef ALLOW_ECCO
16     # include "ECCO_OPTIONS.h"
17     #endif
18 jmc 1.1
19 jmc 1.134 #ifdef ALLOW_AUTODIFF
20 heimbach 1.145 # ifdef ALLOW_GGL90
21     # include "GGL90_OPTIONS.h"
22     # endif
23 jmc 1.1 # ifdef ALLOW_GMREDI
24     # include "GMREDI_OPTIONS.h"
25     # endif
26     # ifdef ALLOW_KPP
27     # include "KPP_OPTIONS.h"
28     # endif
29 jmc 1.29 # ifdef ALLOW_SEAICE
30     # include "SEAICE_OPTIONS.h"
31     # endif
32 gforget 1.117 # ifdef ALLOW_EXF
33     # include "EXF_OPTIONS.h"
34     # endif
35 jmc 1.134 #endif /* ALLOW_AUTODIFF */
36 jmc 1.1
37     CBOP
38     C !ROUTINE: DO_OCEANIC_PHYS
39     C !INTERFACE:
40     SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
41     C !DESCRIPTION: \bv
42     C *==========================================================*
43 jmc 1.28 C | SUBROUTINE DO_OCEANIC_PHYS
44     C | o Controlling routine for oceanic physics and
45 jmc 1.1 C | parameterization
46     C *==========================================================*
47     C | o originally, part of S/R thermodynamics
48     C *==========================================================*
49     C \ev
50    
51 jmc 1.143 C !CALLING SEQUENCE:
52     C DO_OCEANIC_PHYS
53     C |
54     C |-- OBCS_CALC
55     C |
56 jmc 1.147 C |-- OCN_APPLY_IMPORT
57     C |
58 jmc 1.143 C |-- FRAZIL_CALC_RHS
59     C |
60     C |-- THSICE_MAIN
61     C |
62     C |-- SEAICE_FAKE
63     C |-- SEAICE_MODEL
64     C |-- SEAICE_COST_SENSI
65     C |
66 jmc 1.147 C |-- OCN_EXPORT_DATA
67     C |
68 jmc 1.143 C |-- SHELFICE_THERMODYNAMICS
69     C |
70     C |-- ICEFRONT_THERMODYNAMICS
71     C |
72     C |-- SALT_PLUME_DO_EXCH
73     C |
74     C |-- FREEZE_SURFACE
75     C |
76     C |-- EXTERNAL_FORCING_SURF
77     C |
78     C |- k loop (Nr:1):
79     C | - DWNSLP_CALC_RHO
80     C | - BBL_CALC_RHO
81     C | - FIND_RHO_2D @ p(k)
82     C | - FIND_RHO_2D @ p(k-1)
83     C | - GRAD_SIGMA
84     C | - CALC_IVDC
85     C | - DIAGS_RHO_L
86     C |- end k loop.
87     C |
88     C |-- CALC_OCE_MXLAYER
89     C |
90     C |-- SALT_PLUME_CALC_DEPTH
91     C |-- SALT_PLUME_VOLFRAC
92     C |-- SALT_PLUME_APPLY
93     C |-- SALT_PLUME_APPLY
94     C |-- SALT_PLUME_FORCING_SURF
95     C |
96     C |-- KPP_CALC
97     C |-- KPP_CALC_DUMMY
98     C |
99     C |-- PP81_CALC
100     C |
101     C |-- KL10_CALC
102     C |
103     C |-- MY82_CALC
104     C |
105     C |-- GGL90_CALC
106     C |
107     C |-- TIMEAVE_SURF_FLUX
108     C |
109     C |-- GMREDI_CALC_TENSOR
110     C |-- GMREDI_CALC_TENSOR_DUMMY
111     C |
112     C |-- DWNSLP_CALC_FLOW
113     C |-- DWNSLP_CALC_FLOW
114     C |
115 jmc 1.146 C |-- OFFLINE_GET_DIFFUS
116     C |
117 jmc 1.143 C |-- BBL_CALC_RHS
118     C |
119     C |-- MYPACKAGE_CALC_RHS
120     C |
121     C |-- GMREDI_DO_EXCH
122     C |
123     C |-- KPP_DO_EXCH
124     C |
125     C |-- DIAGS_RHO_G
126     C |-- DIAGS_OCEANIC_SURF_FLUX
127     C |-- SALT_PLUME_DIAGNOSTICS_FILL
128     C |
129     C |-- ECCO_PHYS
130    
131 jmc 1.1 C !USES:
132     IMPLICIT NONE
133     C == Global variables ===
134     #include "SIZE.h"
135     #include "EEPARAMS.h"
136     #include "PARAMS.h"
137 jmc 1.69 #include "GRID.h"
138 jmc 1.1 #include "DYNVARS.h"
139 jmc 1.20 #ifdef ALLOW_TIMEAVE
140 jmc 1.131 # include "TIMEAVE_STATV.h"
141     #endif
142     #ifdef ALLOW_OFFLINE
143     # include "OFFLINE_SWITCH.h"
144 jmc 1.20 #endif
145 jmc 1.1
146 jmc 1.134 #ifdef ALLOW_AUTODIFF
147 heimbach 1.80 # include "AUTODIFF_MYFIELDS.h"
148 jmc 1.1 # include "tamc.h"
149     # include "tamc_keys.h"
150     # include "FFIELDS.h"
151 heimbach 1.54 # include "SURFACE.h"
152 jmc 1.1 # include "EOS.h"
153 heimbach 1.145 # ifdef ALLOW_GMREDI
154     # include "GMREDI.h"
155     # endif
156 jmc 1.1 # ifdef ALLOW_KPP
157     # include "KPP.h"
158     # endif
159 gforget 1.91 # ifdef ALLOW_GGL90
160     # include "GGL90.h"
161     # endif
162 jmc 1.1 # ifdef ALLOW_EBM
163     # include "EBM.h"
164     # endif
165 jmc 1.29 # ifdef ALLOW_EXF
166     # include "ctrl.h"
167 jmc 1.40 # include "EXF_FIELDS.h"
168 jmc 1.29 # ifdef ALLOW_BULKFORMULAE
169 jmc 1.40 # include "EXF_CONSTANTS.h"
170 jmc 1.29 # endif
171     # endif
172     # ifdef ALLOW_SEAICE
173 jmc 1.114 # include "SEAICE_SIZE.h"
174 jmc 1.29 # include "SEAICE.h"
175 gforget 1.124 # include "SEAICE_PARAMS.h"
176 jmc 1.29 # endif
177 heimbach 1.105 # ifdef ALLOW_THSICE
178     # include "THSICE_VARS.h"
179     # endif
180 heimbach 1.75 # ifdef ALLOW_SALT_PLUME
181     # include "SALT_PLUME.h"
182     # endif
183 atn 1.142 # ifdef ALLOW_ECCO
184     # ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
185     # include "ecco_cost.h"
186     # endif
187     # endif
188 jmc 1.134 #endif /* ALLOW_AUTODIFF */
189 jmc 1.1
190     C !INPUT/OUTPUT PARAMETERS:
191     C == Routine arguments ==
192 jmc 1.18 C myTime :: Current time in simulation
193     C myIter :: Current iteration number in simulation
194     C myThid :: Thread number for this instance of the routine.
195 jmc 1.1 _RL myTime
196     INTEGER myIter
197     INTEGER myThid
198    
199     C !LOCAL VARIABLES:
200     C == Local variables
201 jmc 1.47 C rhoK, rhoKm1 :: Density at current level, and level above
202 jmc 1.18 C iMin, iMax :: Ranges and sub-block indices on which calculations
203 jmc 1.1 C jMin, jMax are applied.
204 jmc 1.18 C bi, bj :: tile indices
205 jmc 1.147 C msgBuf :: Temp. for building output string
206 jmc 1.18 C i,j,k :: loop indices
207 jmc 1.47 _RL rhoKp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
208     _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
209 jmc 1.1 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
210     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
211     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
212     INTEGER iMin, iMax
213     INTEGER jMin, jMax
214     INTEGER bi, bj
215 jmc 1.18 INTEGER i, j, k
216 jmc 1.147 CHARACTER*(MAX_LEN_MBUF) msgBuf
217 jmc 1.17 INTEGER doDiagsRho
218 jmc 1.131 LOGICAL calcGMRedi, calcKPP, calcConvect
219 jmc 1.17 #ifdef ALLOW_DIAGNOSTICS
220     LOGICAL DIAGNOSTICS_IS_ON
221     EXTERNAL DIAGNOSTICS_IS_ON
222     #endif /* ALLOW_DIAGNOSTICS */
223 jmc 1.141 #ifdef ALLOW_AUTODIFF
224     _RL thetaRef
225     #endif /* ALLOW_AUTODIFF */
226 jmc 1.1 CEOP
227    
228 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
229     C-- dummy statement to end declaration part
230     itdkey = 1
231     #endif /* ALLOW_AUTODIFF_TAMC */
232    
233 jmc 1.1 #ifdef ALLOW_DEBUG
234 jmc 1.96 IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
235 jmc 1.1 #endif
236 jmc 1.36
237 jmc 1.17 doDiagsRho = 0
238     #ifdef ALLOW_DIAGNOSTICS
239     IF ( useDiagnostics .AND. fluidIsWater ) THEN
240 jmc 1.110 IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
241 jmc 1.71 & doDiagsRho = doDiagsRho + 1
242     IF ( DIAGNOSTICS_IS_ON('DRHODR ',myThid) )
243     & doDiagsRho = doDiagsRho + 2
244 jmc 1.110 IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) )
245 jmc 1.71 & doDiagsRho = doDiagsRho + 4
246 jmc 1.110 IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) )
247     & doDiagsRho = doDiagsRho + 8
248 jmc 1.17 ENDIF
249     #endif /* ALLOW_DIAGNOSTICS */
250    
251 jmc 1.131 calcGMRedi = useGMRedi
252     calcKPP = useKPP
253     calcConvect = ivdc_kappa.NE.0.
254     #ifdef ALLOW_OFFLINE
255     IF ( useOffLine ) THEN
256     calcGMRedi = useGMRedi .AND. .NOT.offlineLoadGMRedi
257     calcKPP = useKPP .AND. .NOT.offlineLoadKPP
258     calcConvect=calcConvect.AND. .NOT.offlineLoadConvec
259     ENDIF
260     #endif /* ALLOW_OFFLINE */
261    
262 jmc 1.82 #ifdef ALLOW_OBCS
263     IF (useOBCS) THEN
264     C-- Calculate future values on open boundaries
265     C-- moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
266 heimbach 1.100 # ifdef ALLOW_AUTODIFF_TAMC
267     CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
268     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
269     # endif
270     # ifdef ALLOW_DEBUG
271 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
272 heimbach 1.100 # endif
273 jmc 1.120 CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
274 jahn 1.83 I uVel, vVel, wVel, theta, salt, myThid )
275 jmc 1.82 ENDIF
276     #endif /* ALLOW_OBCS */
277 jmc 1.69
278 jmc 1.147 #ifdef ALLOW_OCN_COMPON_INTERF
279     C-- Apply imported data (from coupled interface) to forcing fields
280     C jmc: moved here before any freezing/seaice pkg adjustment of surf-fluxes
281     IF ( useCoupler ) THEN
282     CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
283     ENDIF
284     #endif /* ALLOW_OCN_COMPON_INTERF */
285    
286 jmc 1.134 #ifdef ALLOW_AUTODIFF
287 gforget 1.87 # ifdef ALLOW_SALT_PLUME
288     DO bj=myByLo(myThid),myByHi(myThid)
289     DO bi=myBxLo(myThid),myBxHi(myThid)
290     DO j=1-OLy,sNy+OLy
291     DO i=1-OLx,sNx+OLx
292     saltPlumeDepth(i,j,bi,bj) = 0. _d 0
293     saltPlumeFlux(i,j,bi,bj) = 0. _d 0
294     ENDDO
295     ENDDO
296     ENDDO
297     ENDDO
298     # endif
299 atn 1.142 # ifdef ALLOW_ECCO
300     # ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
301     DO bj=myByLo(myThid),myByHi(myThid)
302     DO bi=myBxLo(myThid),myBxHi(myThid)
303     DO k=1,Nr
304     DO j=1-OLy,sNy+OLy
305     DO i=1-OLx,sNx+OLx
306     sigmaRfield(i,j,k,bi,bj) = 0. _d 0
307     ENDDO
308     ENDDO
309     ENDDO
310     ENDDO
311     ENDDO
312     # endif
313     # endif
314 jmc 1.134 #endif /* ALLOW_AUTODIFF */
315 gforget 1.87
316 dimitri 1.113 #ifdef ALLOW_FRAZIL
317     IF ( useFRAZIL ) THEN
318     C-- Freeze water in the ocean interior and let it rise to the surface
319 heimbach 1.148 CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
320     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
321 dimitri 1.113 CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
322     ENDIF
323     #endif /* ALLOW_FRAZIL */
324    
325 jmc 1.121 #ifndef OLD_THSICE_CALL_SEQUENCE
326     #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
327     IF ( useThSIce .AND. fluidIsWater ) THEN
328 jmc 1.122 # ifdef ALLOW_AUTODIFF_TAMC
329     CADJ STORE uice,vice = comlev1, key = ikey_dynamics,
330     CADJ & kind = isbyte
331     CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
332     CADJ & kind = isbyte
333     CADJ STORE snowHeight, Tsrf = comlev1, key = ikey_dynamics,
334     CADJ & kind = isbyte
335     CADJ STORE Qice1, Qice2 = comlev1, key = ikey_dynamics,
336     CADJ & kind = isbyte
337     CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics,
338     CADJ & kind = isbyte
339 heimbach 1.125 CADJ STORE hocemxl = comlev1, key = ikey_dynamics,
340     CADJ & kind = isbyte
341 jmc 1.128 CADJ STORE icflxsw = comlev1, key = ikey_dynamics,
342 heimbach 1.126 CADJ & kind = isbyte
343 jmc 1.122 CADJ STORE salt,theta = comlev1, key = ikey_dynamics,
344     CADJ & kind = isbyte
345     CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
346     CADJ & kind = isbyte
347     CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics,
348     CADJ & kind = isbyte
349     CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
350     CADJ & kind = isbyte
351     CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
352     CADJ & kind = isbyte
353     # ifdef NONLIN_FRSURF
354     CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics,
355     CADJ & kind = isbyte
356     # endif
357 heimbach 1.126 # endif /* ALLOW_AUTODIFF_TAMC */
358 jmc 1.121 # ifdef ALLOW_DEBUG
359     IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
360     # endif
361     C-- Step forward Therm.Sea-Ice variables
362     C and modify forcing terms including effects from ice
363     CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
364     CALL THSICE_MAIN( myTime, myIter, myThid )
365     CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
366     ENDIF
367     #endif /* ALLOW_THSICE */
368     #endif /* ndef OLD_THSICE_CALL_SEQUENCE */
369    
370 jmc 1.29 #ifdef ALLOW_SEAICE
371 jmc 1.134 # ifdef ALLOW_AUTODIFF
372 gforget 1.124 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
373     CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte
374     CADJ STORE qnet = comlev1, key=ikey_dynamics, kind=isbyte
375     CADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte
376     CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
377     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
378     #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
379     CADJ STORE evap = comlev1, key=ikey_dynamics, kind=isbyte
380     #endif
381 jmc 1.134 IF ( .NOT.useSEAICE .AND. SEAICEadjMODE .EQ. -1 ) THEN
382 gforget 1.124 CALL SEAICE_FAKE( myTime, myIter, myThid )
383     ENDIF
384     CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
385     CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte
386     CADJ STORE qnet = comlev1, key=ikey_dynamics, kind=isbyte
387     CADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte
388     CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
389     CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
390     #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
391     CADJ STORE evap = comlev1, key=ikey_dynamics, kind=isbyte
392     #endif
393 jmc 1.134 # endif /* ALLOW_AUTODIFF */
394 gforget 1.124 #endif /* ALLOW_SEAICE */
395    
396     #ifdef ALLOW_SEAICE
397 jmc 1.29 IF ( useSEAICE ) THEN
398 heimbach 1.62 # ifdef ALLOW_AUTODIFF_TAMC
399 heimbach 1.65 cph-adj-test(
400 heimbach 1.81 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
401     CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
402 heimbach 1.88 CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte
403 jmc 1.132 CADJ STORE tices = comlev1, key=ikey_dynamics, kind=isbyte
404     CADJ STORE empmr, qnet = comlev1, key=ikey_dynamics, kind=isbyte
405     CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte
406     CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte
407     cCADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
408     cCADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte
409 heimbach 1.65 cph-adj-test)
410 jmc 1.132 c#ifdef ALLOW_EXF
411 heimbach 1.77 CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
412     CADJ & kind = isbyte
413     CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
414     CADJ & kind = isbyte
415     CADJ STORE evap = comlev1, key = ikey_dynamics,
416     CADJ & kind = isbyte
417 jmc 1.132 CADJ STORE uwind,vwind = comlev1, key = ikey_dynamics,
418     CADJ & kind = isbyte
419     c#endif
420 heimbach 1.77 CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
421     CADJ & kind = isbyte
422 heimbach 1.95 # ifdef SEAICE_CGRID
423 heimbach 1.93 CADJ STORE stressdivergencex = comlev1, key = ikey_dynamics,
424     CADJ & kind = isbyte
425     CADJ STORE stressdivergencey = comlev1, key = ikey_dynamics,
426     CADJ & kind = isbyte
427     # endif
428 heimbach 1.62 # ifdef SEAICE_ALLOW_DYNAMICS
429 heimbach 1.77 CADJ STORE uice = comlev1, key = ikey_dynamics,
430     CADJ & kind = isbyte
431     CADJ STORE vice = comlev1, key = ikey_dynamics,
432     CADJ & kind = isbyte
433 jmc 1.132 CADJ STORE dwatn = comlev1, key = ikey_dynamics,
434     CADJ & kind = isbyte
435 heimbach 1.62 # ifdef SEAICE_ALLOW_EVP
436 heimbach 1.77 CADJ STORE seaice_sigma1 = comlev1, key = ikey_dynamics,
437     CADJ & kind = isbyte
438     CADJ STORE seaice_sigma2 = comlev1, key = ikey_dynamics,
439     CADJ & kind = isbyte
440     CADJ STORE seaice_sigma12 = comlev1, key = ikey_dynamics,
441     CADJ & kind = isbyte
442 heimbach 1.62 # endif
443     # endif
444 jmc 1.132 # ifdef SEAICE_VARIABLE_SALINITY
445     CADJ STORE hsalt = comlev1, key = ikey_dynamics,
446 heimbach 1.77 CADJ & kind = isbyte
447 jmc 1.132 # endif
448 heimbach 1.62 # ifdef ATMOSPHERIC_LOADING
449 heimbach 1.77 CADJ STORE pload = comlev1, key = ikey_dynamics,
450     CADJ & kind = isbyte
451     CADJ STORE siceload = comlev1, key = ikey_dynamics,
452     CADJ & kind = isbyte
453 heimbach 1.62 # endif
454     # ifdef NONLIN_FRSURF
455 heimbach 1.77 CADJ STORE recip_hfacc = comlev1, key = ikey_dynamics,
456     CADJ & kind = isbyte
457 heimbach 1.62 # endif
458 heimbach 1.78 # ifdef ANNUAL_BALANCE
459     CADJ STORE balance_itcount = comlev1, key = ikey_dynamics,
460     CADJ & kind = isbyte
461     # endif /* ANNUAL_BALANCE */
462 heimbach 1.126 # ifdef ALLOW_THSICE
463 jmc 1.132 C-- store thSIce vars before advection (called from SEAICE_MODEL) update them:
464     CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
465     CADJ & kind = isbyte
466     CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics,
467 heimbach 1.126 CADJ & kind = isbyte
468 jmc 1.132 CADJ STORE Qice1, Qice2 = comlev1, key = ikey_dynamics,
469     CADJ & kind = isbyte
470     # endif /* ALLOW_THSICE */
471 heimbach 1.126 # endif /* ALLOW_AUTODIFF_TAMC */
472 heimbach 1.62 # ifdef ALLOW_DEBUG
473 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
474 heimbach 1.62 # endif
475 jmc 1.29 CALL TIMER_START('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
476     CALL SEAICE_MODEL( myTime, myIter, myThid )
477     CALL TIMER_STOP ('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
478 heimbach 1.62 # ifdef ALLOW_COST
479 heimbach 1.57 CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
480 heimbach 1.62 # endif
481 heimbach 1.35 ENDIF
482 jmc 1.29 #endif /* ALLOW_SEAICE */
483    
484 jmc 1.147 #if (defined ALLOW_OCN_COMPON_INTERF) && (defined ALLOW_THSICE)
485     C-- After seaice-dyn and advection of pkg/thsice fields,
486     C Export ocean coupling fields to coupled interface (only with pkg/thsice)
487     IF ( useCoupler ) THEN
488     # ifdef ALLOW_DEBUG
489     IF (debugMode) CALL DEBUG_CALL('OCN_EXPORT_DATA',myThid)
490     # endif
491     CALL TIMER_START('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid)
492     CALL OCN_EXPORT_DATA( myTime, myIter, myThid )
493     CALL TIMER_STOP ('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid)
494     ENDIF
495     #endif /* ALLOW_OCN_COMPON_INTERF & ALLOW_THSICE */
496    
497 heimbach 1.64 #ifdef ALLOW_AUTODIFF_TAMC
498 heimbach 1.77 CADJ STORE sst, sss = comlev1, key = ikey_dynamics,
499     CADJ & kind = isbyte
500     CADJ STORE qsw = comlev1, key = ikey_dynamics,
501     CADJ & kind = isbyte
502 heimbach 1.64 # ifdef ALLOW_SEAICE
503 heimbach 1.77 CADJ STORE area = comlev1, key = ikey_dynamics,
504     CADJ & kind = isbyte
505 heimbach 1.64 # endif
506     #endif
507    
508 jmc 1.121 #ifdef OLD_THSICE_CALL_SEQUENCE
509 jscott 1.30 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
510 jmc 1.14 IF ( useThSIce .AND. fluidIsWater ) THEN
511 heimbach 1.101 # ifdef ALLOW_AUTODIFF_TAMC
512     cph(
513     # ifdef NONLIN_FRSURF
514     CADJ STORE uice,vice = comlev1, key = ikey_dynamics,
515     CADJ & kind = isbyte
516     CADJ STORE salt,theta = comlev1, key = ikey_dynamics,
517     CADJ & kind = isbyte
518     CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics,
519     CADJ & kind = isbyte
520     CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics,
521     CADJ & kind = isbyte
522     # endif
523     # endif
524     # ifdef ALLOW_DEBUG
525 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
526 heimbach 1.101 # endif
527 jmc 1.5 C-- Step forward Therm.Sea-Ice variables
528     C and modify forcing terms including effects from ice
529     CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
530     CALL THSICE_MAIN( myTime, myIter, myThid )
531     CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
532     ENDIF
533     #endif /* ALLOW_THSICE */
534 jmc 1.121 #endif /* OLD_THSICE_CALL_SEQUENCE */
535 jmc 1.5
536 mlosch 1.21 #ifdef ALLOW_SHELFICE
537     IF ( useShelfIce .AND. fluidIsWater ) THEN
538     #ifdef ALLOW_DEBUG
539 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
540 mlosch 1.21 #endif
541 heimbach 1.148 #ifdef ALLOW_AUTODIFF_TAMC
542     CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
543     CADJ & kind = isbyte
544     CADJ STORE uvel, vvel = comlev1, key = ikey_dynamics,
545     CADJ & kind = isbyte
546     #endif
547 jmc 1.47 C compute temperature and (virtual) salt flux at the
548 mlosch 1.21 C shelf-ice ocean interface
549     CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
550     & myThid)
551     CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
552     CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
553     & myThid)
554     ENDIF
555     #endif /* ALLOW_SHELFICE */
556    
557 dimitri 1.85 #ifdef ALLOW_ICEFRONT
558     IF ( useICEFRONT .AND. fluidIsWater ) THEN
559     #ifdef ALLOW_DEBUG
560 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
561 dimitri 1.85 #endif
562     C compute temperature and (virtual) salt flux at the
563     C ice-front ocean interface
564     CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
565     & myThid)
566     CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
567     CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
568     & myThid)
569     ENDIF
570     #endif /* ALLOW_ICEFRONT */
571    
572 jmc 1.112 #ifdef ALLOW_SALT_PLUME
573     IF ( useSALT_PLUME ) THEN
574 heimbach 1.136 Catn: exchanging saltPlumeFlux:
575 jmc 1.112 CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
576     ENDIF
577     #endif /* ALLOW_SALT_PLUME */
578    
579 jmc 1.5 C-- Freeze water at the surface
580 heimbach 1.104 IF ( allowFreezing ) THEN
581 jmc 1.5 #ifdef ALLOW_AUTODIFF_TAMC
582 heimbach 1.77 CADJ STORE theta = comlev1, key = ikey_dynamics,
583     CADJ & kind = isbyte
584 jmc 1.5 #endif
585 jmc 1.129 CALL FREEZE_SURFACE( myTime, myIter, myThid )
586 jmc 1.5 ENDIF
587    
588 jmc 1.130 iMin = 1-OLx
589     iMax = sNx+OLx
590     jMin = 1-OLy
591     jMax = sNy+OLy
592    
593     C--- Determines forcing terms based on external fields
594     C relaxation terms, etc.
595 jmc 1.134 #ifdef ALLOW_AUTODIFF
596 jmc 1.130 CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
597     CADJ & kind = isbyte
598 jmc 1.134 #else /* ALLOW_AUTODIFF */
599 jmc 1.130 C-- if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
600     C and all vertical mixing schemes, but keep OBCS_CALC
601 jmc 1.127 IF ( fluidIsWater ) THEN
602 jmc 1.134 #endif /* ALLOW_AUTODIFF */
603 jmc 1.146 #ifdef ALLOW_DEBUG
604     IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
605     #endif
606 jmc 1.130 CALL EXTERNAL_FORCING_SURF(
607     I iMin, iMax, jMin, jMax,
608     I myTime, myIter, myThid )
609 jmc 1.25
610 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
611     C-- HPF directive to help TAMC
612     CHPF$ INDEPENDENT
613     #endif /* ALLOW_AUTODIFF_TAMC */
614     DO bj=myByLo(myThid),myByHi(myThid)
615     #ifdef ALLOW_AUTODIFF_TAMC
616 heimbach 1.15 C-- HPF directive to help TAMC
617     CHPF$ INDEPENDENT
618 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
619     DO bi=myBxLo(myThid),myBxHi(myThid)
620    
621     #ifdef ALLOW_AUTODIFF_TAMC
622     act1 = bi - myBxLo(myThid)
623     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
624     act2 = bj - myByLo(myThid)
625     max2 = myByHi(myThid) - myByLo(myThid) + 1
626     act3 = myThid - 1
627     max3 = nTx*nTy
628     act4 = ikey_dynamics - 1
629     itdkey = (act1 + 1) + act2*max1
630     & + act3*max1*max2
631     & + act4*max1*max2*max3
632 jmc 1.74 #endif /* ALLOW_AUTODIFF_TAMC */
633 jmc 1.1
634     C-- Set up work arrays with valid (i.e. not NaN) values
635     C These inital values do not alter the numerical results. They
636     C just ensure that all memory references are to valid floating
637     C point numbers. This prevents spurious hardware signals due to
638     C uninitialised but inert locations.
639 jmc 1.131 DO k=1,Nr
640     DO j=1-OLy,sNy+OLy
641     DO i=1-OLx,sNx+OLx
642     C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
643     sigmaX(i,j,k) = 0. _d 0
644     sigmaY(i,j,k) = 0. _d 0
645     sigmaR(i,j,k) = 0. _d 0
646     ENDDO
647     ENDDO
648     ENDDO
649 jmc 1.1
650 jmc 1.134 #ifdef ALLOW_AUTODIFF
651 jmc 1.1 DO j=1-OLy,sNy+OLy
652     DO i=1-OLx,sNx+OLx
653 jmc 1.69 rhoKm1 (i,j) = 0. _d 0
654 jmc 1.47 rhoKp1 (i,j) = 0. _d 0
655 jmc 1.1 ENDDO
656     ENDDO
657 jmc 1.131 cph all the following init. are necessary for TAF
658     cph although some of these are re-initialised later.
659 jmc 1.1 DO k=1,Nr
660     DO j=1-OLy,sNy+OLy
661     DO i=1-OLx,sNx+OLx
662 heimbach 1.109 rhoInSitu(i,j,k,bi,bj) = 0.
663 jmc 1.131 # ifdef ALLOW_GGL90
664     GGL90viscArU(i,j,k,bi,bj) = 0. _d 0
665     GGL90viscArV(i,j,k,bi,bj) = 0. _d 0
666     GGL90diffKr(i,j,k,bi,bj) = 0. _d 0
667     # endif /* ALLOW_GGL90 */
668 heimbach 1.136 # ifdef ALLOW_SALT_PLUME
669     # ifdef SALT_PLUME_VOLUME
670     SPforcingS(i,j,k,bi,bj) = 0. _d 0
671     SPforcingT(i,j,k,bi,bj) = 0. _d 0
672     # endif
673     # endif /* ALLOW_SALT_PLUME */
674 jmc 1.131 ENDDO
675     ENDDO
676     ENDDO
677     #ifdef ALLOW_OFFLINE
678     IF ( calcConvect ) THEN
679     #endif
680     DO k=1,Nr
681     DO j=1-OLy,sNy+OLy
682     DO i=1-OLx,sNx+OLx
683 jmc 1.1 IVDConvCount(i,j,k,bi,bj) = 0.
684 jmc 1.131 ENDDO
685     ENDDO
686     ENDDO
687     #ifdef ALLOW_OFFLINE
688     ENDIF
689     IF ( calcGMRedi ) THEN
690     #endif
691 jmc 1.1 # ifdef ALLOW_GMREDI
692 jmc 1.131 DO k=1,Nr
693     DO j=1-OLy,sNy+OLy
694     DO i=1-OLx,sNx+OLx
695 jmc 1.1 Kwx(i,j,k,bi,bj) = 0. _d 0
696     Kwy(i,j,k,bi,bj) = 0. _d 0
697     Kwz(i,j,k,bi,bj) = 0. _d 0
698     # ifdef GM_NON_UNITY_DIAGONAL
699     Kux(i,j,k,bi,bj) = 0. _d 0
700     Kvy(i,j,k,bi,bj) = 0. _d 0
701     # endif
702     # ifdef GM_EXTRA_DIAGONAL
703     Kuz(i,j,k,bi,bj) = 0. _d 0
704     Kvz(i,j,k,bi,bj) = 0. _d 0
705     # endif
706     # ifdef GM_BOLUS_ADVEC
707     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
708     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
709     # endif
710     # ifdef GM_VISBECK_VARIABLE_K
711     VisbeckK(i,j,bi,bj) = 0. _d 0
712     # endif
713 jmc 1.131 ENDDO
714     ENDDO
715     ENDDO
716 jmc 1.1 # endif /* ALLOW_GMREDI */
717 jmc 1.131 #ifdef ALLOW_OFFLINE
718     ENDIF
719     IF ( calcKPP ) THEN
720     #endif
721 heimbach 1.42 # ifdef ALLOW_KPP
722 jmc 1.131 DO k=1,Nr
723     DO j=1-OLy,sNy+OLy
724     DO i=1-OLx,sNx+OLx
725 heimbach 1.42 KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
726     KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
727 jmc 1.1 ENDDO
728     ENDDO
729     ENDDO
730 jmc 1.131 # endif /* ALLOW_KPP */
731     #ifdef ALLOW_OFFLINE
732     ENDIF
733     #endif
734 jmc 1.134 #endif /* ALLOW_AUTODIFF */
735 jmc 1.1
736     #ifdef ALLOW_AUTODIFF_TAMC
737 jmc 1.96 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
738 heimbach 1.77 CADJ & kind = isbyte
739 jmc 1.96 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
740 heimbach 1.77 CADJ & kind = isbyte
741 heimbach 1.12 CADJ STORE totphihyd(:,:,:,bi,bj)
742 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
743 heimbach 1.77 CADJ & kind = isbyte
744 heimbach 1.10 # ifdef ALLOW_KPP
745 jmc 1.96 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
746 heimbach 1.77 CADJ & kind = isbyte
747 jmc 1.96 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
748 heimbach 1.77 CADJ & kind = isbyte
749 heimbach 1.10 # endif
750 heimbach 1.115 # ifdef ALLOW_SALT_PLUME
751     CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
752     CADJ & kind = isbyte
753 heimbach 1.135 CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj, key=itdkey,
754     CADJ & kind = isbyte
755 heimbach 1.115 # endif
756 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
757    
758 jmc 1.71 C-- Always compute density (stored in common block) here; even when it is not
759     C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
760     #ifdef ALLOW_DEBUG
761 jmc 1.110 IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
762 jmc 1.71 #endif
763 jmc 1.134 #ifdef ALLOW_AUTODIFF
764 jmc 1.110 IF ( fluidIsWater ) THEN
765 jmc 1.134 #endif /* ALLOW_AUTODIFF */
766 jmc 1.69 #ifdef ALLOW_DOWN_SLOPE
767 jmc 1.110 IF ( useDOWN_SLOPE ) THEN
768     DO k=1,Nr
769 jmc 1.69 CALL DWNSLP_CALC_RHO(
770     I theta, salt,
771 jmc 1.71 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
772 jmc 1.69 I k, bi, bj, myTime, myIter, myThid )
773 jmc 1.110 ENDDO
774     ENDIF
775 jmc 1.71 #endif /* ALLOW_DOWN_SLOPE */
776 dimitri 1.107 #ifdef ALLOW_BBL
777 jmc 1.110 IF ( useBBL ) THEN
778 dimitri 1.108 C pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
779     C To reduce computation and storage requirement, these densities are stored in the
780     C dry grid boxes of rhoInSitu. See BBL_CALC_RHO for details.
781 jmc 1.110 DO k=Nr,1,-1
782 dimitri 1.107 CALL BBL_CALC_RHO(
783     I theta, salt,
784     O rhoInSitu,
785     I k, bi, bj, myTime, myIter, myThid )
786    
787 jmc 1.110 ENDDO
788     ENDIF
789 dimitri 1.107 #endif /* ALLOW_BBL */
790 jmc 1.110 IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
791     DO k=1,Nr
792 jmc 1.71 CALL FIND_RHO_2D(
793     I iMin, iMax, jMin, jMax, k,
794     I theta(1-OLx,1-OLy,k,bi,bj),
795     I salt (1-OLx,1-OLy,k,bi,bj),
796     O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
797     I k, bi, bj, myThid )
798 jmc 1.110 ENDDO
799     ENDIF
800 jmc 1.134 #ifdef ALLOW_AUTODIFF
801 jmc 1.110 ELSE
802 jmc 1.74 C- fluid is not water:
803 jmc 1.110 DO k=1,Nr
804 jmc 1.141 IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
805     C- isothermal (theta=const) reference state
806     thetaRef = thetaConst
807     ELSE
808     C- horizontally uniform (tRef) reference state
809     thetaRef = tRef(k)
810     ENDIF
811 jmc 1.74 DO j=1-OLy,sNy+OLy
812     DO i=1-OLx,sNx+OLx
813 jmc 1.141 rhoInSitu(i,j,k,bi,bj) =
814     & ( theta(i,j,k,bi,bj)
815     & *( salt(i,j,k,bi,bj)*atm_Rq + oneRL )
816     & - thetaRef )*maskC(i,j,k,bi,bj)
817 jmc 1.74 ENDDO
818     ENDDO
819 jmc 1.110 ENDDO
820     ENDIF
821 jmc 1.134 #endif /* ALLOW_AUTODIFF */
822 jmc 1.110
823     #ifdef ALLOW_DEBUG
824 jmc 1.147 IF (debugMode) THEN
825     WRITE(msgBuf,'(A,2(I4,A))')
826     & 'ENTERING UPWARD K LOOP (bi=', bi, ', bj=', bj,')'
827     CALL DEBUG_MSG(msgBuf(1:43),myThid)
828     ENDIF
829 jmc 1.110 #endif
830    
831     C-- Start of diagnostic loop
832     DO k=Nr,1,-1
833    
834     #ifdef ALLOW_AUTODIFF_TAMC
835     C? Patrick, is this formula correct now that we change the loop range?
836     C? Do we still need this?
837     cph kkey formula corrected.
838     cph Needed for rhoK, rhoKm1, in the case useGMREDI.
839     kkey = (itdkey-1)*Nr + k
840 jmc 1.74 #endif /* ALLOW_AUTODIFF_TAMC */
841 jmc 1.69
842 jmc 1.110 c#ifdef ALLOW_AUTODIFF_TAMC
843     cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
844     cCADJ & kind = isbyte
845     cCADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
846     cCADJ & kind = isbyte
847     c#endif /* ALLOW_AUTODIFF_TAMC */
848    
849 jmc 1.1 C-- Calculate gradients of potential density for isoneutral
850     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
851 jmc 1.131 IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
852 jmc 1.140 & .OR. usePP81 .OR. useKL10
853     & .OR. useMY82 .OR. useGGL90
854 dimitri 1.61 & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
855 jmc 1.1 IF (k.GT.1) THEN
856     #ifdef ALLOW_AUTODIFF_TAMC
857 jmc 1.96 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
858 heimbach 1.77 CADJ & kind = isbyte
859 jmc 1.96 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
860 heimbach 1.77 CADJ & kind = isbyte
861 jmc 1.96 CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey,
862 heimbach 1.77 CADJ & kind = isbyte
863 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
864 jmc 1.68 CALL FIND_RHO_2D(
865     I iMin, iMax, jMin, jMax, k,
866     I theta(1-OLx,1-OLy,k-1,bi,bj),
867     I salt (1-OLx,1-OLy,k-1,bi,bj),
868     O rhoKm1,
869     I k-1, bi, bj, myThid )
870 jmc 1.1 ENDIF
871     #ifdef ALLOW_DEBUG
872 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
873 jmc 1.1 #endif
874 heimbach 1.31 cph Avoid variable aliasing for adjoint !!!
875     DO j=jMin,jMax
876     DO i=iMin,iMax
877 jmc 1.71 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
878 heimbach 1.31 ENDDO
879     ENDDO
880 jmc 1.1 CALL GRAD_SIGMA(
881     I bi, bj, iMin, iMax, jMin, jMax, k,
882 jmc 1.71 I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
883 jmc 1.1 O sigmaX, sigmaY, sigmaR,
884     I myThid )
885 atn 1.142 #ifdef ALLOW_ECCO
886     # ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
887     DO j=jMin,jMax
888     DO i=iMin,iMax
889     sigmaRfield(i,j,k,bi,bj)=sigmaR(i,j,k)
890     ENDDO
891     ENDDO
892     # endif
893     #endif /* ALLOW_ECCO */
894 jmc 1.134 #ifdef ALLOW_AUTODIFF
895 jmc 1.69 #ifdef GMREDI_WITH_STABLE_ADJOINT
896 gforget 1.66 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
897     cgf -> cuts adjoint dependency from slope to state
898 jmc 1.69 CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
899     CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
900     CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
901 gforget 1.66 #endif
902 jmc 1.134 #endif /* ALLOW_AUTODIFF */
903 jmc 1.1 ENDIF
904    
905     C-- Implicit Vertical Diffusion for Convection
906 jmc 1.131 IF (k.GT.1 .AND. calcConvect) THEN
907 jmc 1.1 #ifdef ALLOW_DEBUG
908 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
909 jmc 1.1 #endif
910     CALL CALC_IVDC(
911     I bi, bj, iMin, iMax, jMin, jMax, k,
912 mlosch 1.111 I sigmaR,
913 jmc 1.1 I myTime, myIter, myThid)
914     ENDIF
915    
916 jmc 1.17 #ifdef ALLOW_DIAGNOSTICS
917 jmc 1.110 IF ( doDiagsRho.GE.4 ) THEN
918     CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
919     I rhoInSitu(1-OLx,1-OLy,1,bi,bj),
920 jmc 1.74 I rhoKm1, wVel,
921 jmc 1.71 I myTime, myIter, myThid )
922 jmc 1.17 ENDIF
923     #endif
924    
925 jmc 1.1 C-- end of diagnostic k loop (Nr:1)
926     ENDDO
927    
928 heimbach 1.57 #ifdef ALLOW_AUTODIFF_TAMC
929 jmc 1.69 CADJ STORE IVDConvCount(:,:,:,bi,bj)
930 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
931 heimbach 1.77 CADJ & kind = isbyte
932 heimbach 1.57 #endif
933    
934 jmc 1.47 C-- Diagnose Mixed Layer Depth:
935 jmc 1.131 IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
936 jmc 1.71 CALL CALC_OCE_MXLAYER(
937     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
938     I bi, bj, myTime, myIter, myThid )
939 jmc 1.47 ENDIF
940 heimbach 1.53
941 dimitri 1.52 #ifdef ALLOW_SALT_PLUME
942 dimitri 1.61 IF ( useSALT_PLUME ) THEN
943 jmc 1.71 CALL SALT_PLUME_CALC_DEPTH(
944     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
945     I bi, bj, myTime, myIter, myThid )
946 heimbach 1.136 #ifdef SALT_PLUME_VOLUME
947     CALL SALT_PLUME_VOLFRAC(
948     I bi, bj, myTime, myIter, myThid )
949     C-- get forcings for kpp
950     CALL SALT_PLUME_APPLY(
951     I 1, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
952     I theta, 0,
953     I myTime, myIter, myThid )
954     CALL SALT_PLUME_APPLY(
955     I 2, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
956     I salt, 0,
957     I myTime, myIter, myThid )
958     C-- need to call this S/R from here to apply just before kpp
959     CALL SALT_PLUME_FORCING_SURF(
960     I bi, bj, iMin, iMax, jMin, jMax,
961     I myTime, myIter, myThid )
962     #endif /* SALT_PLUME_VOLUME */
963 dimitri 1.60 ENDIF
964 dimitri 1.61 #endif /* ALLOW_SALT_PLUME */
965    
966 jmc 1.8 #ifdef ALLOW_DIAGNOSTICS
967 jmc 1.71 IF ( MOD(doDiagsRho,4).GE.2 ) THEN
968 jmc 1.16 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
969     & 2, bi, bj, myThid)
970 jmc 1.8 ENDIF
971 dimitri 1.61 #endif /* ALLOW_DIAGNOSTICS */
972 jmc 1.8
973 jmc 1.130 C-- This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called;
974     C now called earlier, before bi,bj loop.
975 jmc 1.1
976     #ifdef ALLOW_AUTODIFF_TAMC
977     cph needed for KPP
978 jmc 1.4 CADJ STORE surfaceForcingU(:,:,bi,bj)
979 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
980 heimbach 1.77 CADJ & kind = isbyte
981 jmc 1.4 CADJ STORE surfaceForcingV(:,:,bi,bj)
982 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
983 heimbach 1.77 CADJ & kind = isbyte
984 jmc 1.4 CADJ STORE surfaceForcingS(:,:,bi,bj)
985 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
986 heimbach 1.77 CADJ & kind = isbyte
987 jmc 1.4 CADJ STORE surfaceForcingT(:,:,bi,bj)
988 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
989 heimbach 1.77 CADJ & kind = isbyte
990 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
991 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
992 heimbach 1.77 CADJ & kind = isbyte
993 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
994    
995     #ifdef ALLOW_KPP
996     C-- Compute KPP mixing coefficients
997 jmc 1.131 IF ( calcKPP ) THEN
998 jmc 1.1 #ifdef ALLOW_DEBUG
999 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
1000 jmc 1.1 #endif
1001 dfer 1.76 CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
1002 jmc 1.1 CALL KPP_CALC(
1003 jmc 1.44 I bi, bj, myTime, myIter, myThid )
1004 dfer 1.76 CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
1005 jmc 1.134 #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
1006 jmc 1.1 ELSE
1007     CALL KPP_CALC_DUMMY(
1008 jmc 1.44 I bi, bj, myTime, myIter, myThid )
1009 jmc 1.134 #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
1010 jmc 1.1 ENDIF
1011     #endif /* ALLOW_KPP */
1012    
1013 mlosch 1.6 #ifdef ALLOW_PP81
1014     C-- Compute PP81 mixing coefficients
1015     IF (usePP81) THEN
1016     #ifdef ALLOW_DEBUG
1017 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
1018 mlosch 1.6 #endif
1019     CALL PP81_CALC(
1020 jmc 1.116 I bi, bj, sigmaR, myTime, myIter, myThid )
1021 mlosch 1.6 ENDIF
1022     #endif /* ALLOW_PP81 */
1023    
1024 jmc 1.140 #ifdef ALLOW_KL10
1025     C-- Compute KL10 mixing coefficients
1026     IF (useKL10) THEN
1027     #ifdef ALLOW_DEBUG
1028     IF (debugMode) CALL DEBUG_CALL('KL10_CALC',myThid)
1029     #endif
1030     CALL KL10_CALC(
1031     I bi, bj, sigmaR, myTime, myIter, myThid )
1032     ENDIF
1033     #endif /* ALLOW_KL10 */
1034    
1035 mlosch 1.6 #ifdef ALLOW_MY82
1036     C-- Compute MY82 mixing coefficients
1037     IF (useMY82) THEN
1038     #ifdef ALLOW_DEBUG
1039 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
1040 mlosch 1.6 #endif
1041     CALL MY82_CALC(
1042 jmc 1.116 I bi, bj, sigmaR, myTime, myIter, myThid )
1043 mlosch 1.6 ENDIF
1044     #endif /* ALLOW_MY82 */
1045    
1046 mlosch 1.9 #ifdef ALLOW_GGL90
1047 gforget 1.91 #ifdef ALLOW_AUTODIFF_TAMC
1048     CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
1049     CADJ & kind = isbyte
1050     #endif /* ALLOW_AUTODIFF_TAMC */
1051 mlosch 1.9 C-- Compute GGL90 mixing coefficients
1052     IF (useGGL90) THEN
1053     #ifdef ALLOW_DEBUG
1054 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
1055 mlosch 1.9 #endif
1056 dfer 1.76 CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
1057 mlosch 1.9 CALL GGL90_CALC(
1058 jmc 1.116 I bi, bj, sigmaR, myTime, myIter, myThid )
1059 dfer 1.76 CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
1060 mlosch 1.9 ENDIF
1061     #endif /* ALLOW_GGL90 */
1062    
1063 jmc 1.20 #ifdef ALLOW_TIMEAVE
1064 jmc 1.36 IF ( taveFreq.GT. 0. _d 0 ) THEN
1065 jmc 1.20 CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
1066     ENDIF
1067 jmc 1.131 IF ( taveFreq.GT.0. .AND. calcConvect ) THEN
1068 jmc 1.20 CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
1069 jmc 1.120 I Nr, deltaTClock, bi, bj, myThid)
1070 jmc 1.20 ENDIF
1071     #endif /* ALLOW_TIMEAVE */
1072    
1073 jmc 1.69 #ifdef ALLOW_GMREDI
1074 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
1075     # ifndef GM_EXCLUDE_CLIPPING
1076     cph storing here is needed only for one GMREDI_OPTIONS:
1077     cph define GM_BOLUS_ADVEC
1078     cph keep it although TAF says you dont need to.
1079 jmc 1.86 cph but I have avoided the #ifdef for now, in case more things change
1080 jmc 1.96 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey,
1081 heimbach 1.77 CADJ & kind = isbyte
1082 jmc 1.96 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey,
1083 heimbach 1.77 CADJ & kind = isbyte
1084 jmc 1.96 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey,
1085 heimbach 1.77 CADJ & kind = isbyte
1086 jmc 1.47 # endif
1087     #endif /* ALLOW_AUTODIFF_TAMC */
1088    
1089     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
1090 jmc 1.131 IF ( calcGMRedi ) THEN
1091 jmc 1.47 #ifdef ALLOW_DEBUG
1092 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
1093 jmc 1.47 #endif
1094     CALL GMREDI_CALC_TENSOR(
1095 jmc 1.51 I iMin, iMax, jMin, jMax,
1096 jmc 1.47 I sigmaX, sigmaY, sigmaR,
1097 jmc 1.51 I bi, bj, myTime, myIter, myThid )
1098 jmc 1.134 #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
1099 jmc 1.47 ELSE
1100     CALL GMREDI_CALC_TENSOR_DUMMY(
1101 jmc 1.51 I iMin, iMax, jMin, jMax,
1102 jmc 1.47 I sigmaX, sigmaY, sigmaR,
1103 jmc 1.51 I bi, bj, myTime, myIter, myThid )
1104 jmc 1.134 #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
1105 jmc 1.47 ENDIF
1106 jmc 1.69 #endif /* ALLOW_GMREDI */
1107    
1108     #ifdef ALLOW_DOWN_SLOPE
1109     IF ( useDOWN_SLOPE ) THEN
1110     C-- Calculate Downsloping Flow for Down_Slope parameterization
1111     IF ( usingPCoords ) THEN
1112     CALL DWNSLP_CALC_FLOW(
1113 jmc 1.71 I bi, bj, kSurfC, rhoInSitu,
1114 jmc 1.69 I myTime, myIter, myThid )
1115     ELSE
1116     CALL DWNSLP_CALC_FLOW(
1117 jmc 1.71 I bi, bj, kLowC, rhoInSitu,
1118 jmc 1.69 I myTime, myIter, myThid )
1119     ENDIF
1120     ENDIF
1121     #endif /* ALLOW_DOWN_SLOPE */
1122 jmc 1.47
1123 jmc 1.106 C-- end bi,bj loops.
1124     ENDDO
1125     ENDDO
1126    
1127 jmc 1.134 #ifndef ALLOW_AUTODIFF
1128 jmc 1.98 C--- if fluid Is Water: end
1129 jmc 1.106 ENDIF
1130 jmc 1.98 #endif
1131    
1132 jmc 1.146 #ifdef ALLOW_OFFLINE
1133     IF ( useOffLine ) THEN
1134     #ifdef ALLOW_DEBUG
1135     IF (debugMode) CALL DEBUG_CALL('OFFLINE_GET_DIFFUS',myThid)
1136     #endif /* ALLOW_DEBUG */
1137     CALL OFFLINE_GET_DIFFUS( myTime, myIter, myThid )
1138     ENDIF
1139     #endif /* ALLOW_OFFLINE */
1140    
1141 dimitri 1.107 #ifdef ALLOW_BBL
1142     IF ( useBBL ) THEN
1143     CALL BBL_CALC_RHS(
1144     I myTime, myIter, myThid )
1145     ENDIF
1146     #endif /* ALLOW_BBL */
1147    
1148 dimitri 1.94 #ifdef ALLOW_MYPACKAGE
1149 jmc 1.106 IF ( useMYPACKAGE ) THEN
1150     CALL MYPACKAGE_CALC_RHS(
1151     I myTime, myIter, myThid )
1152     ENDIF
1153 dimitri 1.94 #endif /* ALLOW_MYPACKAGE */
1154    
1155 jmc 1.99 #ifdef ALLOW_GMREDI
1156 jmc 1.131 IF ( calcGMRedi ) THEN
1157 jmc 1.99 CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
1158     ENDIF
1159     #endif /* ALLOW_GMREDI */
1160    
1161 jmc 1.98 #ifdef ALLOW_KPP
1162 jmc 1.131 IF ( calcKPP ) THEN
1163 jmc 1.45 CALL KPP_DO_EXCH( myThid )
1164     ENDIF
1165 jmc 1.98 #endif /* ALLOW_KPP */
1166 jmc 1.45
1167 jmc 1.18 #ifdef ALLOW_DIAGNOSTICS
1168     IF ( fluidIsWater .AND. useDiagnostics ) THEN
1169 jmc 1.74 CALL DIAGS_RHO_G(
1170 jmc 1.110 I rhoInSitu, uVel, vVel, wVel,
1171 jmc 1.71 I myTime, myIter, myThid )
1172 jmc 1.133 ENDIF
1173     IF ( useDiagnostics ) THEN
1174 jmc 1.18 CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
1175     ENDIF
1176 jmc 1.131 IF ( calcConvect .AND. useDiagnostics ) THEN
1177 jmc 1.71 CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
1178     & 0, Nr, 0, 1, 1, myThid )
1179 jmc 1.19 ENDIF
1180 heimbach 1.136 #ifdef ALLOW_SALT_PLUME
1181     IF ( useDiagnostics )
1182     & CALL SALT_PLUME_DIAGNOSTICS_FILL(bi,bj,myThid)
1183     #endif
1184 jmc 1.18 #endif
1185    
1186 gforget 1.123 #ifdef ALLOW_ECCO
1187 gforget 1.144 IF ( useECCO ) CALL ECCO_PHYS( myThid )
1188 gforget 1.123 #endif
1189    
1190 jmc 1.1 #ifdef ALLOW_DEBUG
1191 jmc 1.96 IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
1192 jmc 1.1 #endif
1193    
1194     RETURN
1195     END

  ViewVC Help
Powered by ViewVC 1.1.22