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