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

Diff of /MITgcm/model/src/thermodynamics.F

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

revision 1.10 by adcroft, Thu Sep 27 18:06:43 2001 UTC revision 1.36 by heimbach, Tue Jan 21 19:19:45 2003 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF_TAMC
6    # ifdef ALLOW_GMREDI
7    #  include "GMREDI_OPTIONS.h"
8    # endif
9    # ifdef ALLOW_KPP
10    #  include "KPP_OPTIONS.h"
11    # endif
12    cswdice --- add ----
13    #ifdef ALLOW_THERM_SEAICE
14    #include "ICE.h"
15    #endif
16    cswdice ------
17    #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
20  C     !ROUTINE: THERMODYNAMICS  C     !ROUTINE: THERMODYNAMICS
# Line 75  C     == Global variables === Line 88  C     == Global variables ===
88  # include "tamc.h"  # include "tamc.h"
89  # include "tamc_keys.h"  # include "tamc_keys.h"
90  # include "FFIELDS.h"  # include "FFIELDS.h"
91    # include "EOS.h"
92  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
93  #  include "KPP.h"  #  include "KPP.h"
94  # endif  # endif
# Line 151  C     This is currently used by IVDC and Line 165  C     This is currently used by IVDC and
165        INTEGER i, j        INTEGER i, j
166        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
167    
 Cjmc : add for phiHyd output <- but not working if multi tile per CPU  
 c     CHARACTER*(MAX_LEN_MBUF) suff  
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
168  CEOP  CEOP
169    
170  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
171  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
172        ikey = 1        ikey = 1
173          itdkey = 1
174  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
175    
176  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
# Line 174  C     uninitialised but inert locations. Line 184  C     uninitialised but inert locations.
184          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
185          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
186          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
         rhoKM1 (i,j) = 0. _d 0  
187          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
188          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
189          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
# Line 208  CHPF$&                  ) Line 211  CHPF$&                  )
211  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
212            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
213            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
214            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
215            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
216            act3 = myThid - 1            act3 = myThid - 1
217            max3 = nTx*nTy            max3 = nTx*nTy
   
218            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
219              itdkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
220       &                      + act3*max1*max2       &                      + act3*max1*max2
221       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
222  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 232  C--     Set up work arrays that need val Line 231  C--     Set up work arrays that need val
231            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
232            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
233            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
234              rhoKM1 (i,j)   = 0. _d 0
235           ENDDO           ENDDO
236          ENDDO          ENDDO
237    
# Line 239  C--     Set up work arrays that need val Line 239  C--     Set up work arrays that need val
239           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
240            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
241  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
242               phiHyd(i,j,k) = 0. _d 0
243               sigmaX(i,j,k) = 0. _d 0
244               sigmaY(i,j,k) = 0. _d 0
245               sigmaR(i,j,k) = 0. _d 0
246             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
247             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
248             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
249  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
250             gT(i,j,k,bi,bj) = 0. _d 0  cph all the following init. are necessary for TAF
251             gS(i,j,k,bi,bj) = 0. _d 0  cph although some of these are re-initialised later.
252  #ifdef ALLOW_PASSIVE_TRACER             gT(i,j,k,bi,bj)   = 0. _d 0
253               gS(i,j,k,bi,bj)   = 0. _d 0
254    # ifdef ALLOW_PASSIVE_TRACER
255             gTr1(i,j,k,bi,bj) = 0. _d 0             gTr1(i,j,k,bi,bj) = 0. _d 0
256  #endif  # endif
257  #endif  # ifdef ALLOW_GMREDI
258               Kwx(i,j,k,bi,bj)  = 0. _d 0
259               Kwy(i,j,k,bi,bj)  = 0. _d 0
260               Kwz(i,j,k,bi,bj)  = 0. _d 0
261    #  ifdef GM_NON_UNITY_DIAGONAL
262               Kux(i,j,k,bi,bj)  = 0. _d 0
263               Kvy(i,j,k,bi,bj)  = 0. _d 0
264    #  endif
265    #  ifdef GM_EXTRA_DIAGONAL
266               Kuz(i,j,k,bi,bj)  = 0. _d 0
267               Kvz(i,j,k,bi,bj)  = 0. _d 0
268    #  endif
269    #  ifdef GM_BOLUS_ADVEC
270               GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
271               GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
272    #  endif
273    #  ifdef GM_VISBECK_VARIABLE_K
274               VisbeckK(i,j,bi,bj)   = 0. _d 0
275    #  endif
276    # endif /* ALLOW_GMREDI */
277    #endif /* ALLOW_AUTODIFF_TAMC */
278            ENDDO            ENDDO
279           ENDDO           ENDDO
280          ENDDO          ENDDO
281    
282          iMin = 1-OLx+1          iMin = 1-OLx
283          iMax = sNx+OLx          iMax = sNx+OLx
284          jMin = 1-OLy+1          jMin = 1-OLy
285          jMax = sNy+OLy          jMax = sNy+OLy
286    
287    
288  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
289  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
290  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
291    #ifdef ALLOW_KPP
292    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
293    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
294    #endif
295  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
296    
297  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 272  C? Patrick, is this formula correct now Line 302  C? Patrick, is this formula correct now
302  C? Do we still need this?  C? Do we still need this?
303  cph kkey formula corrected.  cph kkey formula corrected.
304  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
305           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
306  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
307    
308  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
309            CALL INTEGRATE_FOR_W(  c         CALL INTEGRATE_FOR_W(
310       I                         bi, bj, k, uVel, vVel,  c    I                         bi, bj, k, uVel, vVel,
311       O                         wVel,  c    O                         wVel,
312       I                         myThid )  c    I                         myThid )
313    
314  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
315  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
316  C--       Apply OBC to W if in N-H mode  C--       Apply OBC to W if in N-H mode
317            IF (useOBCS.AND.nonHydrostatic) THEN  c         IF (useOBCS.AND.nonHydrostatic) THEN
318              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
319            ENDIF  c         ENDIF
320  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
321  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
322    
323    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
324    C--       MOST of THERMODYNAMICS will be disabled
325    #ifndef SINGLE_LAYER_MODE
326    
327  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
328  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
329  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
# Line 299  c         IF ( k.GT.1 .AND. (useGMRedi.O Line 331  c         IF ( k.GT.1 .AND. (useGMRedi.O
331  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
332  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
333  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
334    CADJ STORE pressure(:,:,k,bi,bj) =
335    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
336  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
337              CALL FIND_RHO(              CALL FIND_RHO(
338       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
339       I        theta, salt,       I        theta, salt,
340       O        rhoK,       O        rhoK,
341       I        myThid )       I        myThid )
342    
343              IF (k.GT.1) THEN              IF (k.GT.1) THEN
344  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
345  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
346  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
347    CADJ STORE pressure(:,:,k-1,bi,bj) =
348    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
349  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
350               CALL FIND_RHO(               CALL FIND_RHO(
351       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
352       I        theta, salt,       I        theta, salt,
353       O        rhoKm1,       O        rhoKm1,
354       I        myThid )       I        myThid )
# Line 323  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 360  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
360       I             myThid )       I             myThid )
361            ENDIF            ENDIF
362    
363    #ifdef ALLOW_AUTODIFF_TAMC
364    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
365    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
366    #endif /* ALLOW_AUTODIFF_TAMC */
367  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
368  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
369            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
# Line 333  c ==> should use sigmaR !!! Line 374  c ==> should use sigmaR !!!
374       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
375            ENDIF            ENDIF
376    
377    #endif /* SINGLE_LAYER_MODE */
378    
379  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
380          ENDDO          ENDDO
381    
382  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
383  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
384  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
385    CADJ STORE pressure (:,:,:,bi,bj) =
386    CADJ &     comlev1_bibj, key=itdkey, byte=isbyte
387  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
388    
389  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
390  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
391          IF (useOBCS) THEN          IF (useOBCS) THEN
392            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
393       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
394       I            myThid )       I            myThid )
395          ENDIF          ENDIF
396  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
397    
398    
399    c********************************************
400    cswdice --- add ---
401    #ifdef ALLOW_THERM_SEAICE
402    C--     Determines forcing terms based on external fields
403    c--     including effects from ice
404            CALL ICE_FORCING(
405         I             bi, bj, iMin, iMax, jMin, jMax,
406         I             myThid )
407    #else
408    
409    cswdice --- end add ---
410    
411  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
412  C       relaxation terms, etc.  C       relaxation terms, etc.
413          CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
414       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
415       I             myThid )       I             myThid )
416    cswdice --- add ----
417    #endif
418    cswdice --- end add ---
419    c******************************************
420    
421    
422    
423    
424  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
425  cph needed for KPP  cph needed for KPP
426  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
427  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
428  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
429  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
430  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
431  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
432  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
433  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
434  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
435    
436    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
437    C--     MOST of THERMODYNAMICS will be disabled
438    #ifndef SINGLE_LAYER_MODE
439    
440  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
441    
442  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
443  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
444  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  cph define GM_BOLUS_ADVEC
445  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
446    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
447    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
448    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
449  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
450    
451  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
452          IF (useGMRedi) THEN          IF (useGMRedi) THEN
453            DO k=1,Nr            CALL GMREDI_CALC_TENSOR(
454              CALL GMREDI_CALC_TENSOR(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
455       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
456       I             myThid )       I             myThid )
           ENDDO  
457  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
458          ELSE          ELSE
459            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
460              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
461       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
462       I             myThid )       I             myThid )
           ENDDO  
463  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
464          ENDIF          ENDIF
465    
466  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
467  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
468  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
469  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
470  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
471    
472  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 415  C--     Compute KPP mixing coefficients Line 485  C--     Compute KPP mixing coefficients
485    
486  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
487  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE KPPghat   (:,:,:,bi,bj)
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
488  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
489  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
490  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
491  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
492  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
493    
494  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
495    
496  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
497  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
498  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
499  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
500  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
501  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
502  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
503  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
504  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
505  #endif  #endif
506  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
507    
508  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
509  C       AIM - atmospheric intermediate model, physics package code.  C       AIM - atmospheric intermediate model, physics package code.
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
510          IF ( useAIM ) THEN          IF ( useAIM ) THEN
511           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
512           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
513           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
514          ENDIF          ENDIF
515  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
516    
517    #ifdef ALLOW_TIMEAVE
518            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
519              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
520         I                               deltaTclock, bi, bj, myThid)
521            ENDIF
522    #endif /* ALLOW_TIMEAVE */
523    
524    #ifndef DISABLE_MULTIDIM_ADVECTION
525  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
526  C       method in the absence of any other terms and, if used, is done here.  C       method in the absence of any other terms and, if used, is done here.
527          IF (multiDimAdvection) THEN  C
528           IF (tempStepping .AND.  C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
529       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.  C The default is to use multi-dimensinal advection for non-linear advection
530       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.  C schemes. However, for the sake of efficiency of the adjoint it is necessary
531       &       tempAdvScheme.NE.ENUM_CENTERED_4TH )  C to be able to exclude this scheme to avoid excessive storage and
532       &   CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,theta,  C recomputation. It *is* differentiable, if you need it.
533       U                      gT,  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
534    C disable this section of code.
535            IF (tempMultiDimAdvec) THEN
536              CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
537         U                      theta,gT,
538       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
539           IF (saltStepping .AND.          ENDIF
540       &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.          IF (saltMultiDimAdvec) THEN
541       &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
542       &       saltAdvScheme.NE.ENUM_CENTERED_4TH )       U                      salt,gS,
      &   CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,salt,  
      U                      gS,  
543       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
544          ENDIF          ENDIF
545    C Since passive tracers are configurable separately from T,S we
546    C call the multi-dimensional method for PTRACERS regardless
547    C of whether multiDimAdvection is set or not.
548    #ifdef ALLOW_PTRACERS
549            IF ( usePTRACERS ) THEN
550             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
551            ENDIF
552    #endif /* ALLOW_PTRACERS */
553    #endif /* DISABLE_MULTIDIM_ADVECTION */
554    
555  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
556          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 472  C--     Start of thermodynamics loop Line 558  C--     Start of thermodynamics loop
558  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
559  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
560  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
561           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
562  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
563    
564  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 494  C--       Get temporary terms used by te Line 580  C--       Get temporary terms used by te
580       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
581       I         myThid)       I         myThid)
582    
583    #ifdef ALLOW_GMREDI
584    
585    C--   Residual transp = Bolus transp + Eulerian transp
586              IF (useGMRedi) THEN
587                CALL GMREDI_CALC_UVFLOW(
588         &            uTrans, vTrans, bi, bj, k, myThid)
589                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
590         &                    rTrans, bi, bj, k, myThid)
591              ENDIF
592    
593    #ifdef ALLOW_AUTODIFF_TAMC
594    #ifdef GM_BOLUS_ADVEC
595    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
596    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
597    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
598    #endif
599    #endif /* ALLOW_AUTODIFF_TAMC */
600    
601    #endif /* ALLOW_GMREDI */
602    
603  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
604  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
605  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 525  C        and step forward storing result Line 631  C        and step forward storing result
631             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
632       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
633       I         theta, gT,       I         theta, gT,
      U         gTnm1,  
634       I         myIter, myThid)       I         myIter, myThid)
635           ENDIF           ENDIF
636    cswdice ---- add ---
637    #ifdef ALLOW_THERM_SEAICE
638           if (k.eq.1) then
639            call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )
640           endif
641    #endif
642    cswdice -- end add ---
643           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
644             CALL CALC_GS(             CALL CALC_GS(
645       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
# Line 538  C        and step forward storing result Line 650  C        and step forward storing result
650             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
651       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
652       I         salt, gS,       I         salt, gS,
      U         gSnm1,  
653       I         myIter, myThid)       I         myIter, myThid)
654           ENDIF           ENDIF
655  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
# Line 552  C        and step forward storing result Line 663  C        and step forward storing result
663             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
664       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
665       I         Tr1, gTr1,       I         Tr1, gTr1,
      U         gTr1NM1,  
666       I         myIter,myThid)       I         myIter,myThid)
667           ENDIF           ENDIF
668  #endif  #endif
669    #ifdef ALLOW_PTRACERS
670             IF ( usePTRACERS ) THEN
671               CALL PTRACERS_INTEGERATE(
672         I         bi,bj,k,
673         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
674         X         KappaRS,
675         I         myIter,myTime,myThid)
676             ENDIF
677    #endif /* ALLOW_PTRACERS */
678    
679  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
680  C--      Apply open boundary conditions  C--      Apply open boundary conditions
# Line 565  C--      Apply open boundary conditions Line 684  C--      Apply open boundary conditions
684  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
685    
686  C--      Freeze water  C--      Freeze water
687           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN
688  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
689  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
690  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
691  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
692              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
# Line 576  CADJ &   , key = kkey, byte = isbyte Line 695  CADJ &   , key = kkey, byte = isbyte
695  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
696          ENDDO          ENDDO
697    
698    cswdice -- add ---
699    #ifdef ALLOW_THERM_SEAICE
700    c timeaveraging for ice model values
701               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
702    #endif
703    cswdice --- end add ---
704    
705    
706    
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey dont seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isnt needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
707    
708  C--     Implicit diffusion  C--     Implicit diffusion
709          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
710    
711           IF (tempStepping) THEN           IF (tempStepping) THEN
712  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
713              idkey = iikey + 1  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
714  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
715              CALL IMPLDIFF(              CALL IMPLDIFF(
716       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 604  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_ Line 721  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_
721    
722           IF (saltStepping) THEN           IF (saltStepping) THEN
723  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
724           idkey = iikey + 2  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
725  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
726              CALL IMPLDIFF(              CALL IMPLDIFF(
727       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 617  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_ Line 733  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_
733  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
734           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
735  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
736  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
737  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
738            CALL IMPLDIFF(            CALL IMPLDIFF(
739       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
# Line 627  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev Line 743  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev
743           ENDIF           ENDIF
744  #endif  #endif
745    
746    #ifdef ALLOW_PTRACERS
747    C Vertical diffusion (implicit) for passive tracers
748             IF ( usePTRACERS ) THEN
749               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
750             ENDIF
751    #endif /* ALLOW_PTRACERS */
752    
753  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
754  C--      Apply open boundary conditions  C--      Apply open boundary conditions
755           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 639  C--      Apply open boundary conditions Line 762  C--      Apply open boundary conditions
762  C--     End If implicitDiffusion  C--     End If implicitDiffusion
763          ENDIF          ENDIF
764    
765    #endif /* SINGLE_LAYER_MODE */
766    
767  Ccs-  Ccs-
768         ENDDO         ENDDO
769        ENDDO        ENDDO
770    
771  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
772        IF ( useAIM ) THEN  c     IF ( useAIM ) THEN
773         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
774        ENDIF  c     ENDIF
        _EXCH_XYZ_R8(gT,myThid)  
        _EXCH_XYZ_R8(gS,myThid)  
 #else  
       IF (staggerTimeStep.AND.useCubedSphereExchange) THEN  
        _EXCH_XYZ_R8(gT,myThid)  
        _EXCH_XYZ_R8(gS,myThid)  
       ENDIF  
775  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
776    c     IF ( staggerTimeStep ) THEN
777    c      IF ( useAIM .OR. useCubedSphereExchange ) THEN
778    c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
779    c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
780    c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
781    cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
782    c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
783    c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
784    c      ENDIF  
785    c     ENDIF  
786    
787    #ifndef DISABLE_DEBUGMODE
788          If (debugMode) THEN
789           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
790           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
791           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
792           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
793           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
794           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
795           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
796           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
797           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
798    #ifdef ALLOW_PTRACERS
799           IF ( usePTRACERS ) THEN
800             CALL PTRACERS_DEBUG(myThid)
801           ENDIF
802    #endif /* ALLOW_PTRACERS */
803          ENDIF
804    #endif
805    
806        RETURN        RETURN
807        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22