/[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.16 by jmc, Fri Feb 8 22:15:33 2002 UTC revision 1.34 by heimbach, Fri Jan 10 19:06:05 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 156  CEOP Line 170  CEOP
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 169  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 216  CHPF$&                  )
216            act3 = myThid - 1            act3 = myThid - 1
217            max3 = nTx*nTy            max3 = nTx*nTy
218            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
219            ikey = (act1 + 1) + act2*max1            itdkey = (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 223  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 230  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    # endif /* ALLOW_GMREDI */
274    #endif /* ALLOW_AUTODIFF_TAMC */
275            ENDDO            ENDDO
276           ENDDO           ENDDO
277          ENDDO          ENDDO
278    
279          iMin = 1-OLx+1          iMin = 1-OLx
280          iMax = sNx+OLx          iMax = sNx+OLx
281          jMin = 1-OLy+1          jMin = 1-OLy
282          jMax = sNy+OLy          jMax = sNy+OLy
283    
284    
285  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
286  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
287  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
288  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
289  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
290  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
291  #endif  #endif
292  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
293    
# Line 267  C? Patrick, is this formula correct now Line 299  C? Patrick, is this formula correct now
299  C? Do we still need this?  C? Do we still need this?
300  cph kkey formula corrected.  cph kkey formula corrected.
301  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
302           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  
303  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
304    
305  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
306            CALL INTEGRATE_FOR_W(  c         CALL INTEGRATE_FOR_W(
307       I                         bi, bj, k, uVel, vVel,  c    I                         bi, bj, k, uVel, vVel,
308       O                         wVel,  c    O                         wVel,
309       I                         myThid )  c    I                         myThid )
310    
311  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
312  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
313  C--       Apply OBC to W if in N-H mode  C--       Apply OBC to W if in N-H mode
314            IF (useOBCS.AND.nonHydrostatic) THEN  c         IF (useOBCS.AND.nonHydrostatic) THEN
315              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
316            ENDIF  c         ENDIF
317  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
318  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
319    
320    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
321    C--       MOST of THERMODYNAMICS will be disabled
322    #ifndef SINGLE_LAYER_MODE
323    
324  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
325  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
326  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 294  c         IF ( k.GT.1 .AND. (useGMRedi.O Line 328  c         IF ( k.GT.1 .AND. (useGMRedi.O
328  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
329  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
330  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
331    CADJ STORE pressure(:,:,k,bi,bj) =
332    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
333  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
334              CALL FIND_RHO(              CALL FIND_RHO(
335       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
336       I        theta, salt,       I        theta, salt,
337       O        rhoK,       O        rhoK,
338       I        myThid )       I        myThid )
339    
340              IF (k.GT.1) THEN              IF (k.GT.1) THEN
341  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
342  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
343  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
344    CADJ STORE pressure(:,:,k-1,bi,bj) =
345    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
346  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
347               CALL FIND_RHO(               CALL FIND_RHO(
348       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
349       I        theta, salt,       I        theta, salt,
350       O        rhoKm1,       O        rhoKm1,
351       I        myThid )       I        myThid )
# Line 318  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 357  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
357       I             myThid )       I             myThid )
358            ENDIF            ENDIF
359    
360    #ifdef ALLOW_AUTODIFF_TAMC
361    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
362    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
363    #endif /* ALLOW_AUTODIFF_TAMC */
364  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
365  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
366            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
# Line 328  c ==> should use sigmaR !!! Line 371  c ==> should use sigmaR !!!
371       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
372            ENDIF            ENDIF
373    
374    #endif /* SINGLE_LAYER_MODE */
375    
376  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
377          ENDDO          ENDDO
378    
379  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
380  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
381  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
382    CADJ STORE pressure (:,:,:,bi,bj) =
383    CADJ &     comlev1_bibj, key=itdkey, byte=isbyte
384  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
385    
386  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
# Line 345  C--     Calculate future values on open Line 392  C--     Calculate future values on open
392          ENDIF          ENDIF
393  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
394    
395    
396    c********************************************
397    cswdice --- add ---
398    #ifdef ALLOW_THERM_SEAICE
399    C--     Determines forcing terms based on external fields
400    c--     including effects from ice
401            CALL ICE_FORCING(
402         I             bi, bj, iMin, iMax, jMin, jMax,
403         I             myThid )
404    #else
405    
406    cswdice --- end add ---
407    
408  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
409  C       relaxation terms, etc.  C       relaxation terms, etc.
410          CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
411       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
412       I             myThid )       I             myThid )
413    cswdice --- add ----
414    #endif
415    cswdice --- end add ---
416    c******************************************
417    
418    
419    
420    
421  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
422  cph needed for KPP  cph needed for KPP
423  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
424  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
425  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
426  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
427  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
428  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
429  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
430  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
431  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
432    
433    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
434    C--     MOST of THERMODYNAMICS will be disabled
435    #ifndef SINGLE_LAYER_MODE
436    
437  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
438    
439  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
440  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
441  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  cph define GM_BOLUS_ADVEC
442  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
443    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
444    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
445    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
446  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
447  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
448          IF (useGMRedi) THEN          IF (useGMRedi) THEN
# Line 385  C--     Calculate iso-neutral slopes for Line 460  C--     Calculate iso-neutral slopes for
460          ENDIF          ENDIF
461    
462  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
463  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
464  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
465  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
466  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
467    
468  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 409  CADJ STORE KPPghat   (:,:,:,bi,bj) Line 484  CADJ STORE KPPghat   (:,:,:,bi,bj)
484  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
485  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
486  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
487  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
488  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
489    
490  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
491    
492  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
493  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
494  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
495  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
496  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
497  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
498  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
499  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
500  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
501  #endif  #endif
502  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
503    
504  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
505  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  
506          IF ( useAIM ) THEN          IF ( useAIM ) THEN
507           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
508           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
509           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
510          ENDIF          ENDIF
511  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
512    
# Line 454  C to be able to exclude this scheme to a Line 528  C to be able to exclude this scheme to a
528  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
529  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
530  C disable this section of code.  C disable this section of code.
531          IF (multiDimAdvection) THEN          IF (tempMultiDimAdvec) THEN
          IF (tempStepping .AND.  
      &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.  
      &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.  
      &       tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  
532            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
533       U                      theta,gT,       U                      theta,gT,
534       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
535           ENDIF          ENDIF
536           IF (saltStepping .AND.          IF (saltMultiDimAdvec) THEN
      &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.  
      &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.  
      &       saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  
537            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
538       U                      salt,gS,       U                      salt,gS,
539       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
          ENDIF  
540          ENDIF          ENDIF
541    C Since passive tracers are configurable separately from T,S we
542    C call the multi-dimensional method for PTRACERS regardless
543    C of whether multiDimAdvection is set or not.
544    #ifdef ALLOW_PTRACERS
545            IF ( usePTRACERS ) THEN
546             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
547            ENDIF
548    #endif /* ALLOW_PTRACERS */
549  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
550    
551  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
# Line 480  C--     Start of thermodynamics loop Line 554  C--     Start of thermodynamics loop
554  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
555  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
556  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
557           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
558  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
559    
560  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 502  C--       Get temporary terms used by te Line 576  C--       Get temporary terms used by te
576       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
577       I         myThid)       I         myThid)
578    
579    #ifdef ALLOW_GMREDI
580    C--   Residual transp = Bolus transp + Eulerian transp
581              IF (useGMRedi) THEN
582                CALL GMREDI_CALC_UVFLOW(
583         &            uTrans, vTrans, bi, bj, k, myThid)
584                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
585         &                    rTrans, bi, bj, k, myThid)
586              ENDIF
587    #endif /* ALLOW_GMREDI */
588    
589  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
590  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
591  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 535  C        and step forward storing result Line 619  C        and step forward storing result
619       I         theta, gT,       I         theta, gT,
620       I         myIter, myThid)       I         myIter, myThid)
621           ENDIF           ENDIF
622    cswdice ---- add ---
623    #ifdef ALLOW_THERM_SEAICE
624           if (k.eq.1) then
625            call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )
626           endif
627    #endif
628    cswdice -- end add ---
629           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
630             CALL CALC_GS(             CALL CALC_GS(
631       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
# Line 561  C        and step forward storing result Line 652  C        and step forward storing result
652       I         myIter,myThid)       I         myIter,myThid)
653           ENDIF           ENDIF
654  #endif  #endif
655    #ifdef ALLOW_PTRACERS
656             IF ( usePTRACERS ) THEN
657               CALL PTRACERS_INTEGERATE(
658         I         bi,bj,k,
659         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
660         X         KappaRS,
661         I         myIter,myTime,myThid)
662             ENDIF
663    #endif /* ALLOW_PTRACERS */
664    
665  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
666  C--      Apply open boundary conditions  C--      Apply open boundary conditions
# Line 570  C--      Apply open boundary conditions Line 670  C--      Apply open boundary conditions
670  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
671    
672  C--      Freeze water  C--      Freeze water
673           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN
674  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
675  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
676  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
# Line 581  CADJ &   , key = kkey, byte = isbyte Line 681  CADJ &   , key = kkey, byte = isbyte
681  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
682          ENDDO          ENDDO
683    
684    cswdice -- add ---
685    #ifdef ALLOW_THERM_SEAICE
686    c timeaveraging for ice model values
687               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
688    #endif
689    cswdice --- end add ---
690    
691    
692    
 #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 */  
693    
694  C--     Implicit diffusion  C--     Implicit diffusion
695          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
696    
697           IF (tempStepping) THEN           IF (tempStepping) THEN
698  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
699              idkey = iikey + 1  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
700  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
701              CALL IMPLDIFF(              CALL IMPLDIFF(
702       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 609  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib Line 707  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib
707    
708           IF (saltStepping) THEN           IF (saltStepping) THEN
709  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
710           idkey = iikey + 2  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
711  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
712              CALL IMPLDIFF(              CALL IMPLDIFF(
713       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 622  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bib Line 719  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bib
719  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
720           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
721  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
722  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
723  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
724            CALL IMPLDIFF(            CALL IMPLDIFF(
725       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
# Line 632  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b Line 729  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b
729           ENDIF           ENDIF
730  #endif  #endif
731    
732    #ifdef ALLOW_PTRACERS
733    C Vertical diffusion (implicit) for passive tracers
734             IF ( usePTRACERS ) THEN
735               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
736             ENDIF
737    #endif /* ALLOW_PTRACERS */
738    
739  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
740  C--      Apply open boundary conditions  C--      Apply open boundary conditions
741           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 644  C--      Apply open boundary conditions Line 748  C--      Apply open boundary conditions
748  C--     End If implicitDiffusion  C--     End If implicitDiffusion
749          ENDIF          ENDIF
750    
751    #endif /* SINGLE_LAYER_MODE */
752    
753  Ccs-  Ccs-
754         ENDDO         ENDDO
755        ENDDO        ENDDO
756    
757  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
758        IF ( useAIM ) THEN  c     IF ( useAIM ) THEN
759         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
760        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  
761  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
762    c     IF ( staggerTimeStep ) THEN
763    c      IF ( useAIM .OR. useCubedSphereExchange ) THEN
764    c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
765    c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
766    c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
767    cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
768    c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
769    c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
770    c      ENDIF  
771    c     ENDIF  
772    
773    #ifndef DISABLE_DEBUGMODE
774          If (debugMode) THEN
775           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
776           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
777           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
778           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
779           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
780           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
781           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
782           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
783           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
784    #ifdef ALLOW_PTRACERS
785           IF ( usePTRACERS ) THEN
786             CALL PTRACERS_DEBUG(myThid)
787           ENDIF
788    #endif /* ALLOW_PTRACERS */
789          ENDIF
790    #endif
791    
792        RETURN        RETURN
793        END        END
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22