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

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

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

revision 1.11 by heimbach, Mon Sep 27 14:56:42 2004 UTC revision 1.24 by mlosch, Fri Mar 3 10:20:34 2006 UTC
# Line 35  C     == Global variables === Line 35  C     == Global variables ===
35  #include "PARAMS.h"  #include "PARAMS.h"
36  #include "DYNVARS.h"  #include "DYNVARS.h"
37  #include "GRID.h"  #include "GRID.h"
38  c #include "GAD.h"  #ifdef ALLOW_TIMEAVE
39  c #ifdef ALLOW_PTRACERS  #include "TIMEAVE_STATV.h"
40  c #include "PTRACERS_SIZE.h"  #endif
41  c #include "PTRACERS.h"  #if defined (ALLOW_BALANCE_FLUXES) && !(defined ALLOW_AUTODIFF_TAMC)
42  c #endif  #include "FFIELDS.h"
43  c #ifdef ALLOW_TIMEAVE  #endif
 c #include "TIMEAVE_STATV.h"  
 c #endif  
44    
45  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
46  # include "tamc.h"  # include "tamc.h"
# Line 65  c #endif Line 63  c #endif
63    
64  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
65  C     == Routine arguments ==  C     == Routine arguments ==
66  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
67  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
68  C     myThid - Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
69        _RL myTime        _RL myTime
70        INTEGER myIter        INTEGER myIter
71        INTEGER myThid        INTEGER myThid
72    
73  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
74  C     == Local variables  C     == Local variables
75  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKM1  :: Density at current level, and level above
76  C     useVariableK   = T when vertical diffusion is not constant  C     iMin, iMax    :: Ranges and sub-block indices on which calculations
 C     iMin, iMax     - Ranges and sub-block indices on which calculations  
77  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
78  C     bi, bj  C     bi, bj        :: tile indices
79  C     k, kup,        - Index for layer above and below. kup and kDown  C     i,j,k         :: loop indices
 C     kDown, km1       are switched with layer to be the appropriate  
 C                      index into fVerTerm.  
80        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
84        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL kp1Msk  
       LOGICAL useVariableK  
85        INTEGER iMin, iMax        INTEGER iMin, iMax
86        INTEGER jMin, jMax        INTEGER jMin, jMax
87        INTEGER bi, bj        INTEGER bi, bj
88        INTEGER i, j        INTEGER i, j, k
89        INTEGER k, km1, kup, kDown        INTEGER doDiagsRho
90        INTEGER iTracer, ip  #ifdef ALLOW_DIAGNOSTICS
91          LOGICAL  DIAGNOSTICS_IS_ON
92          EXTERNAL DIAGNOSTICS_IS_ON
93    #endif /* ALLOW_DIAGNOSTICS */
94    
95  CEOP  CEOP
96    
97    #ifdef ALLOW_AUTODIFF_TAMC
98    C--   dummy statement to end declaration part
99          itdkey = 1
100    #endif /* ALLOW_AUTODIFF_TAMC */
101    
102  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
103        IF ( debugLevel .GE. debLevB )        IF ( debugLevel .GE. debLevB )
104       &    CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)       &    CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
105  #endif  #endif
106    
107          doDiagsRho = 0
108    #ifdef ALLOW_DIAGNOSTICS
109          IF ( useDiagnostics .AND. fluidIsWater ) THEN
110            IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) ) doDiagsRho = 1
111            IF ( DIAGNOSTICS_IS_ON('RHOANOSQ',myThid) .OR.
112         &       DIAGNOSTICS_IS_ON('URHOMASS',myThid) .OR.
113         &       DIAGNOSTICS_IS_ON('VRHOMASS',myThid) .OR.
114         &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) .OR.
115         &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) ) doDiagsRho = 2
116          ENDIF
117    #endif /* ALLOW_DIAGNOSTICS */
118    
119  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
120        IF ( useThSIce .AND. buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN        IF ( useThSIce .AND. fluidIsWater ) THEN
121  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
122          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
123       &    CALL DEBUG_CALL('THSICE_MAIN',myThid)       &    CALL DEBUG_CALL('THSICE_MAIN',myThid)
# Line 117  C       and modify forcing terms includi Line 130  C       and modify forcing terms includi
130        ENDIF        ENDIF
131  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
132    
133    #ifdef ALLOW_SHELFICE
134          IF ( useShelfIce .AND. fluidIsWater ) THEN
135    #ifdef ALLOW_DEBUG
136            IF ( debugLevel .GE. debLevB )
137         &    CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
138    #endif
139    C     compute temperature and (virtual) salt flux at the
140    C     shelf-ice ocean interface
141           CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
142         &       myThid)
143           CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
144           CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
145         &      myThid)
146          ENDIF
147    #endif /* ALLOW_SHELFICE */
148    
149  C--   Freeze water at the surface  C--   Freeze water at the surface
150  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
151  CADJ STORE theta = comlev1, key = ikey_dynamics  CADJ STORE theta = comlev1, key = ikey_dynamics
152  #endif  #endif
153        IF ( allowFreezing .AND. .NOT. useSEAICE        IF ( allowFreezing
154         &                   .AND. .NOT. useSEAICE
155       &                   .AND. .NOT. useThSIce ) THEN       &                   .AND. .NOT. useThSIce ) THEN
156          CALL FREEZE_SURFACE(  myTime, myIter, myThid )          CALL FREEZE_SURFACE(  myTime, myIter, myThid )
157        ENDIF        ENDIF
# Line 137  C jmc: do not know precisely where to pu Line 167  C jmc: do not know precisely where to pu
167  #endif /* COMPONENT_MODULE */  #endif /* COMPONENT_MODULE */
168    
169  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 C--   dummy statement to end declaration part  
       ikey = 1  
       itdkey = 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
170  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
171  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
172  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
173        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
   
174  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
175  C--    HPF directive to help TAMC  C--   HPF directive to help TAMC
176  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$ INDEPENDENT
 CHPF$&                  ,utrans,vtrans,xA,yA  
 CHPF$&                  )  
177  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
178         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
179    
180  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 228  cph although some of these are re-initia Line 247  cph although some of these are re-initia
247  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
248  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
249  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
250  CADJ STORE totphihyd  CADJ STORE totphihyd(:,:,:,bi,bj)
251  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
252  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
253  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
254  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
255  # endif  # endif
 # ifdef EXACT_CONSERV  
 CADJ STORE pmepr(:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
256  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
257    
258  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
# Line 258  cph Needed for rhok, rhokm1, in the case Line 274  cph Needed for rhok, rhokm1, in the case
274  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
275  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
276  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
277            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
278         &                   .OR. doDiagsRho.GE.1 ) THEN
279  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
280              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
281       &       CALL DEBUG_CALL('FIND_RHO',myThid)       &       CALL DEBUG_CALL('FIND_RHO',myThid)
# Line 296  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 313  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
313            ENDIF            ENDIF
314    
315  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
316  # ifndef GM_EXCLUDE_CLIPPING  ctest# ifndef GM_EXCLUDE_CLIPPING
317  CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
318  # endif  ctest# endif
319  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
320  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
321  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
# Line 314  c ==> should use sigmaR !!! Line 331  c ==> should use sigmaR !!!
331       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
332            ENDIF            ENDIF
333    
334    #ifdef ALLOW_DIAGNOSTICS
335              IF ( doDiagsRho.GE.2 ) THEN
336                CALL DIAGS_RHO( k, bi, bj,
337         I                      rhoK, rhoKm1,
338         I                      myTime, myIter, myThid)
339              ENDIF
340    #endif
341    
342  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
343          ENDDO          ENDDO
344    
345  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
346          IF ( usediagnostics .AND.  c       IF ( useDiagnostics .AND.
347       &       (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN  c    &       (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN
348            CALL fill_diagnostics (myThid, 'DRHODR  ', 0, Nr,          IF ( doDiagsRho.GE.1 ) THEN
349       &         3, bi, bj, sigmaR)            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
350         &         2, bi, bj, myThid)
351          ENDIF          ENDIF
352  #endif  #endif
353    
# Line 339  C--     Calculate future values on open Line 365  C--     Calculate future values on open
365  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
366    
367  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
368          IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN          IF ( fluidIsWater ) THEN
369  #endif  #endif
370    #ifdef ALLOW_BALANCE_FLUXES
371    C     balance fluxes
372             IF ( balanceEmPmR )
373         &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,
374         &        'EmPmR', myTime, myThid )
375             IF ( balanceQnet )
376         &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,
377         &        'Qnet ', myTime, myThid )
378    #endif /* ALLOW_BALANCE_FLUXES */
379  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
380  C       relaxation terms, etc.  C       relaxation terms, etc.
381  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
382          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
383       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
384  #endif  #endif
385    #ifdef ALLOW_AUTODIFF_TAMC
386    CADJ STORE EmPmR(:,:,bi,bj)
387    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
388    CADJ STORE PmEpR(:,:,bi,bj)
389    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
390    #endif
391           CALL EXTERNAL_FORCING_SURF(           CALL EXTERNAL_FORCING_SURF(
392       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
393       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
# Line 364  CADJ STORE surfaceForcingS(:,:,bi,bj) Line 405  CADJ STORE surfaceForcingS(:,:,bi,bj)
405  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
406  CADJ STORE surfaceForcingT(:,:,bi,bj)  CADJ STORE surfaceForcingT(:,:,bi,bj)
407  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
 # ifdef ALLOW_SEAICE  
408  CADJ STORE surfaceForcingTice(:,:,bi,bj)  CADJ STORE surfaceForcingTice(:,:,bi,bj)
409  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
 # endif  
410  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
411    
412  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
# Line 380  cph keep it although TAF says you dont n Line 419  cph keep it although TAF says you dont n
419  cph but I've avoided the #ifdef for now, in case more things change  cph but I've avoided the #ifdef for now, in case more things change
420  CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
421  CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
422  cnewCADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
423  # endif  # endif
424  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
425    
# Line 459  C--     Compute GGL90 mixing coefficient Line 498  C--     Compute GGL90 mixing coefficient
498          ENDIF          ENDIF
499  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
500    
501    #ifdef ALLOW_TIMEAVE
502            IF ( taveFreq.GT. 0. _d 0 .AND. fluidIsWater ) THEN
503              CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
504            ENDIF
505            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
506              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
507         I                           Nr, deltaTclock, bi, bj, myThid)
508            ENDIF
509    #endif /* ALLOW_TIMEAVE */
510    
511  C--   end bi,bj loops.  C--   end bi,bj loops.
512         ENDDO         ENDDO
513        ENDDO        ENDDO
514    
515    #ifdef ALLOW_DIAGNOSTICS
516          IF ( fluidIsWater .AND. useDiagnostics ) THEN
517            CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
518          ENDIF
519          IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
520            CALL DIAGNOSTICS_FILL( IVDConvCount,'CONVADJ ',
521         &                         0, Nr, 0, 1, 1, myThid )
522          ENDIF
523    #endif
524    
525  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
526           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
527       &    CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)       &    CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22