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

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

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

revision 1.49 by heimbach, Fri Jun 9 02:45:04 2000 UTC revision 1.52 by adcroft, Thu Jun 29 18:49:50 2000 UTC
# Line 36  C     == Global variables === Line 36  C     == Global variables ===
36  #include "DYNVARS.h"  #include "DYNVARS.h"
37  #include "GRID.h"  #include "GRID.h"
38    
 #ifdef ALLOW_KPP  
 #include "KPPMIX.h"  
 #endif  
   
39  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
40  #include "tamc.h"  #include "tamc.h"
41  #include "tamc_keys.h"  #include "tamc_keys.h"
# Line 89  C                      surface height Line 85  C                      surface height
85  C                      anomaly.  C                      anomaly.
86  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     etaSurfX,      - Holds surface elevation gradient in X and Y.
87  C     etaSurfY  C     etaSurfY
 C     K13, K23, K33  - Non-zero elements of small-angle approximation  
 C                      diffusion tensor.  
 C     KapGM          - Spatially varying Visbeck et. al mixing coeff.  
88  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
89  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
90  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
# Line 128  C                      index into fVerTe Line 121  C                      index into fVerTe
121        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL K13     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL K23     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL K33     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL KapGM   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
124        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
125        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
126        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
127        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
128          _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
129          _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
130          _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
131    
132  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
133    C #ifdef INCLUDE_CONVECT_CALL
134        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
135  #endif  C #endif
136    
137        INTEGER iMin, iMax        INTEGER iMin, iMax
138        INTEGER jMin, jMax        INTEGER jMin, jMax
# Line 154  C                      index into fVerTe Line 147  C                      index into fVerTe
147    
148        INTEGER act1, act2, act3, act4        INTEGER act1, act2, act3, act4
149        INTEGER max1, max2, max3        INTEGER max1, max2, max3
150        INTEGER ikact, iikey,kkey        INTEGER iikey, kkey
151        INTEGER maximpl        INTEGER maximpl
152  #endif  #endif
153    
# Line 228  C     uninitialised but inert locations. Line 221  C     uninitialised but inert locations.
221          fMer(i,j)    = 0. _d 0          fMer(i,j)    = 0. _d 0
222          DO K=1,Nr          DO K=1,Nr
223           phiHyd (i,j,k)  = 0. _d 0           phiHyd (i,j,k)  = 0. _d 0
          K13(i,j,k)  = 0. _d 0  
          K23(i,j,k)  = 0. _d 0  
          K33(i,j,k)  = 0. _d 0  
224           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
225           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
226             sigmaX(i,j,k) = 0. _d 0
227             sigmaY(i,j,k) = 0. _d 0
228             sigmaR(i,j,k) = 0. _d 0
229          ENDDO          ENDDO
230          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
231          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
# Line 255  C--   HPF directive to help TAMC Line 248  C--   HPF directive to help TAMC
248  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
249  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
250  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV
251  !HPF$&                  ,phiHyd,K13,K23,K33,KapGM  !HPF$&                  ,phiHyd,
252  !HPF$&                  ,utrans,vtrans,maskc,xA,yA  !HPF$&                  ,utrans,vtrans,maskc,xA,yA
253  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
254  !HPF$&                  )  !HPF$&                  )
# Line 295  C--     Set up work arrays that need val Line 288  C--     Set up work arrays that need val
288            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
289            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
290            phiHyd(i,j,1) = 0. _d 0            phiHyd(i,j,1) = 0. _d 0
           K13   (i,j,1) = 0. _d 0  
           K23   (i,j,1) = 0. _d 0  
           K33   (i,j,1) = 0. _d 0  
           KapGM (i,j)   = GMkbackground  
291           ENDDO           ENDDO
292          ENDDO          ENDDO
293    
# Line 377  CADJ STORE salt (:,:,k,bi,bj)  = comlev1 Line 366  CADJ STORE salt (:,:,k,bi,bj)  = comlev1
366  #endif  #endif
367    
368          IF (       (.NOT. BOTTOM_LAYER)          IF (       (.NOT. BOTTOM_LAYER)
 #ifdef ALLOW_KPP  
      &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
369       &     ) THEN       &     ) THEN
370  C--      Check static stability with layer below  C--      Check static stability with layer below
371  C--      and mix as needed.  C--      and mix as needed.
# Line 439  C--     phiHyd(z=0)=0 Line 425  C--     phiHyd(z=0)=0
425       I      bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1,       I      bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1,
426       U      phiHyd,       U      phiHyd,
427       I      myThid )       I      myThid )
428            CALL GRAD_SIGMA(
429         I            bi, bj, iMin, iMax, jMin, jMax, K,
430         I            rhoKm1, rhoKm1, rhoKm1,
431         O            sigmaX, sigmaY, sigmaR,
432         I            myThid )
433    
434  C----------------------------------------------  C--     Start of downward loop
 C--     start of downward loop  
 C----------------------------------------------  
435          DO K=2,Nr          DO K=2,Nr
436    
437  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
438           kkey = ikact*(Nr-2+1) + (k-2) + 1           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
439  #endif  #endif
440    
441           BOTTOM_LAYER = K .EQ. Nr           BOTTOM_LAYER = K .EQ. Nr
# Line 483  CADJ STORE salt (:,:,k,bi,bj)  = comlev1 Line 472  CADJ STORE salt (:,:,k,bi,bj)  = comlev1
472       I      myThid )       I      myThid )
473  #endif  #endif
474           IF (       (.NOT. BOTTOM_LAYER)           IF (       (.NOT. BOTTOM_LAYER)
 #ifdef ALLOW_KPP  
      &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
475       &      ) THEN       &      ) THEN
476  C--       Check static stability with layer below and mix as needed.  C--       Check static stability with layer below and mix as needed.
477  C--       Density of K+1 level (below W(K+1)) reference to K level.  C--       Density of K+1 level (below W(K+1)) reference to K level.
# Line 554  C--      Calculate iso-neutral slopes fo Line 540  C--      Calculate iso-neutral slopes fo
540       O        rhoTmp,       O        rhoTmp,
541       I        myThid )       I        myThid )
542  #endif  #endif
543             CALL GRAD_SIGMA(
544         I             bi, bj, iMin, iMax, jMin, jMax, K,
545         I             rhoK, rhotmp, rhoK,
546         O             sigmaX, sigmaY, sigmaR,
547         I             myThid )
548    
 #ifdef  INCLUDE_CALC_ISOSLOPES_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE kapgm (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
          CALL CALC_ISOSLOPES(  
      I        bi, bj, iMin, iMax, jMin, jMax, K,  
      I        rhoKm1, rhoK, rhotmp,  
      O        K13, K23, K33, KapGM,  
      I        myThid )  
 #endif  
549    
550           DO J=jMin,jMax           DO J=jMin,jMax
551            DO I=iMin,iMax            DO I=iMin,iMax
# Line 580  CADJ STORE kapgm (:,:)  = comlev1_3d, ke Line 558  CADJ STORE kapgm (:,:)  = comlev1_3d, ke
558          ENDDO          ENDDO
559  C--     end of k loop  C--     end of k loop
560    
561    #ifdef ALLOW_GMREDI
562    #ifdef ALLOW_AUTODIFF_TAMC
563    CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte
564    CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte
565    CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte
566    #endif
567            DO K=1, Nr
568             IF (use_GMRedi) CALL GMREDI_CALC_TENSOR(
569         I             bi, bj, iMin, iMax, jMin, jMax, K,
570         I             sigmaX, sigmaY, sigmaR,
571         I             myThid )
572            ENDDO
573    #endif
574    
575  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
576  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte
577  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte
# Line 588  CADJ STORE vvel (:,:,:,bi,bj)  = comlev1 Line 580  CADJ STORE vvel (:,:,:,bi,bj)  = comlev1
580  #endif  #endif
581    
582  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
 C----------------------------------------------  
583  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
584  C----------------------------------------------          CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)
585          IF (usingKPPmixing) THEN          CALL KPP_CALC(
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE fu  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE fv  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
          CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'  
      I          , myThid)  
          CALL KVMIX(  
586       I               bi, bj, myTime, myThid )       I               bi, bj, myTime, myThid )
587           CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'          CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)
      I        , myThid)  
         ENDIF  
588  #endif  #endif
589    
590  C----------------------------------------------  C--     Start of upward loop
 C--     start of upward loop  
 C----------------------------------------------  
591          DO K = Nr, 1, -1          DO K = Nr, 1, -1
592    
593           kM1  =max(1,k-1)   ! Points to level above k (=k-1)           kM1  =max(1,k-1)   ! Points to level above k (=k-1)
# Line 620  C--------------------------------------- Line 600  C---------------------------------------
600           jMax = sNy+OLy-1           jMax = sNy+OLy-1
601    
602  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
603           kkey = ikact*(Nr-1+1) + (k-1) + 1           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
604  #endif  #endif
605    
606  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 646  C--      Get temporary terms used by ten Line 626  C--      Get temporary terms used by ten
626  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
627           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
628       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,K,
629       I        maskC,maskUp,KapGM,K33,       I        maskC,maskUp,
630       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
631       I        myThid)       I        myThid)
632  #endif  #endif
# Line 676  C--      Calculate active tracer tendenc Line 656  C--      Calculate active tracer tendenc
656            CALL CALC_GT(            CALL CALC_GT(
657       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
658       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
659       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
660       U         aTerm,xTerm,fZon,fMer,fVerT,       U         aTerm,xTerm,fZon,fMer,fVerT,
661       I         myTime, myThid)       I         myTime, myThid)
662           ENDIF           ENDIF
# Line 684  C--      Calculate active tracer tendenc Line 664  C--      Calculate active tracer tendenc
664            CALL CALC_GS(            CALL CALC_GS(
665       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
666       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
667       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
668       U         aTerm,xTerm,fZon,fMer,fVerS,       U         aTerm,xTerm,fZon,fMer,fVerS,
669       I         myTime, myThid)       I         myTime, myThid)
670           ENDIF           ENDIF
# Line 731  C--      Cumulative diagnostic calculati Line 711  C--      Cumulative diagnostic calculati
711           IF (taveFreq.GT.0.) THEN           IF (taveFreq.GT.0.) THEN
712            CALL DO_TIME_AVERAGES(            CALL DO_TIME_AVERAGES(
713       I                           myTime, myIter, bi, bj, K, kUp, kDown,       I                           myTime, myIter, bi, bj, K, kUp, kDown,
714       I                           K13, K23, rVel, KapGM, ConvectCount,       I                           rVel, ConvectCount,
715       I                           myThid )       I                           myThid )
716           ENDIF           ENDIF
717  #endif  #endif
# Line 739  C--      Cumulative diagnostic calculati Line 719  C--      Cumulative diagnostic calculati
719    
720          ENDDO ! K          ENDDO ! K
721    
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
   
722  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
723             maximpl = 6             maximpl = 6
724             iikey = ikact*maximpl             iikey = (ikey-1)*maximpl
725  #endif  #endif
726    
727    C--     Implicit diffusion
728            IF (implicitDiffusion) THEN
729    
730           IF (tempStepping) THEN           IF (tempStepping) THEN
731  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
732              idkey = iikey + 1              idkey = iikey + 1
# Line 831  C    &            maxval(rVel(1:sNx,1:sN Line 811  C    &            maxval(rVel(1:sNx,1:sN
811  C     write(0,*) 'dynamics: rVel(2) ',  C     write(0,*) 'dynamics: rVel(2) ',
812  C    &            minval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.),  C    &            minval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.),
813  C    &            maxval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.)  C    &            maxval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.)
 cblk  write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),  
 cblk &                           maxval(K13(1:sNx,1:sNy,:))  
 cblk  write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),  
 cblk &                           maxval(K23(1:sNx,1:sNy,:))  
 cblk  write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),  
 cblk &                           maxval(K33(1:sNx,1:sNy,:))  
814  C     write(0,*) 'dynamics: gT ',minval(gT(1:sNx,1:sNy,:,:,:)),  C     write(0,*) 'dynamics: gT ',minval(gT(1:sNx,1:sNy,:,:,:)),
815  C    &                           maxval(gT(1:sNx,1:sNy,:,:,:))  C    &                           maxval(gT(1:sNx,1:sNy,:,:,:))
816  C     write(0,*) 'dynamics: T  ',minval(Theta(1:sNx,1:sNy,:,:,:)),  C     write(0,*) 'dynamics: T  ',minval(Theta(1:sNx,1:sNy,:,:,:)),
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.52
  ViewVC Help
Powered by ViewVC 1.1.22