/[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.28 by jmc, Thu Mar 30 02:33:05 2006 UTC revision 1.40 by jmc, Mon Apr 16 23:31:59 2007 UTC
# Line 11  C $Name$ Line 11  C $Name$
11  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
12  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
13  # endif  # endif
14    # ifdef ALLOW_SEAICE
15    #  include "SEAICE_OPTIONS.h"
16    # endif
17  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
# Line 59  C     == Global variables === Line 62  C     == Global variables ===
62  # ifdef EXACT_CONSERV  # ifdef EXACT_CONSERV
63  #  include "SURFACE.h"  #  include "SURFACE.h"
64  # endif  # endif
65    # ifdef ALLOW_EXF
66    #  include "ctrl.h"
67    #  include "EXF_FIELDS.h"
68    #  include "EXF_CLIM_FIELDS.h"
69    #  ifdef ALLOW_BULKFORMULAE
70    #   include "EXF_CONSTANTS.h"
71    #  endif
72    # endif
73    # ifdef ALLOW_SEAICE
74    #  include "SEAICE.h"
75    # endif
76  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
77    
78  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 77  C     iMin, iMax    :: Ranges and sub-bl Line 91  C     iMin, iMax    :: Ranges and sub-bl
91  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
92  C     bi, bj        :: tile indices  C     bi, bj        :: tile indices
93  C     i,j,k         :: loop indices  C     i,j,k         :: loop indices
94          _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 100  C--   dummy statement to end declaration Line 115  C--   dummy statement to end declaration
115  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
116    
117  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
118        IF ( debugLevel .GE. debLevB )        IF ( debugLevel .GE. debLevB )
119       &    CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)       &     CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
120  #endif  #endif
121    
122        doDiagsRho = 0        doDiagsRho = 0
123  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
124        IF ( useDiagnostics .AND. fluidIsWater ) THEN        IF ( useDiagnostics .AND. fluidIsWater ) THEN
# Line 116  C--   dummy statement to end declaration Line 131  C--   dummy statement to end declaration
131        ENDIF        ENDIF
132  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
133    
134  #ifdef ALLOW_THSICE  #ifdef ALLOW_SEAICE
135    C--   Call sea ice model to compute forcing/external data fields.  In
136    C     addition to computing prognostic sea-ice variables and diagnosing the
137    C     forcing/external data fields that drive the ocean model, SEAICE_MODEL
138    C     also sets theta to the freezing point under sea-ice.  The implied
139    C     surface heat flux is then stored in variable surfaceTendencyTice,
140    C     which is needed by KPP package (kpp_calc.F and kpp_transport_t.F)
141    C     to diagnose surface buoyancy fluxes and for the non-local transport
142    C     term.  Because this call precedes model thermodynamics, temperature
143    C     under sea-ice may not be "exactly" at the freezing point by the time
144    C     theta is dumped or time-averaged.
145          IF ( useSEAICE ) THEN
146    #ifdef ALLOW_AUTODIFF_TAMC
147    CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics
148    CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics
149    CADJ STORE theta               = comlev1, key = ikey_dynamics
150    cph# ifdef EXF_READ_EVAP
151    CADJ STORE evap                = comlev1, key = ikey_dynamics
152    cph# endif
153    cph# ifdef SEAICE_ALLOW_DYNAMICS
154    CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics
155    cph# endif
156    # ifdef SEAICE_CGRID
157    CADJ STORE fu, fv              = comlev1, key = ikey_dynamics
158    CADJ STORE seaicemasku         = comlev1, key = ikey_dynamics
159    CADJ STORE seaicemaskv         = comlev1, key = ikey_dynamics
160    #  ifdef ATMOSPHERIC_LOADING
161    CADJ STORE siceload            = comlev1, key = ikey_dynamics
162    #  endif
163    # endif
164    #endif
165    #ifdef ALLOW_DEBUG
166            IF ( debugLevel .GE. debLevB )
167         &    CALL DEBUG_CALL('SEAICE_MODEL',myThid)
168    #endif
169            CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
170            CALL SEAICE_MODEL( myTime, myIter, myThid )
171            CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
172    #ifdef ALLOW_COST_ICE
173            CALL COST_ICE_TEST ( myTime, myIter, myThid )
174    #endif
175          ENDIF
176    #endif /* ALLOW_SEAICE */
177    
178    #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
179        IF ( useThSIce .AND. fluidIsWater ) THEN        IF ( useThSIce .AND. fluidIsWater ) THEN
180  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
181          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
# Line 159  CADJ STORE theta = comlev1, key = ikey_d Line 218  CADJ STORE theta = comlev1, key = ikey_d
218  #ifdef ALLOW_OCN_COMPON_INTERF  #ifdef ALLOW_OCN_COMPON_INTERF
219  C--    Apply imported data (from coupled interface) to forcing fields  C--    Apply imported data (from coupled interface) to forcing fields
220  C jmc: do not know precisely where to put this call (bf or af thSIce ?)  C jmc: do not know precisely where to put this call (bf or af thSIce ?)
221         IF ( useCoupler ) THEN        IF ( useCoupler ) THEN
222           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
223         ENDIF        ENDIF
224  #endif /* ALLOW_OCN_COMPON_INTERF */  #endif /* ALLOW_OCN_COMPON_INTERF */
225    
226  #ifdef ALLOW_BALANCE_FLUXES  #ifdef ALLOW_BALANCE_FLUXES
227  C     balance fluxes  C     balance fluxes
228         IF ( balanceEmPmR )        IF ( balanceEmPmR )
229       &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,       &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,
230       &        'EmPmR', myTime, myThid )       &        'EmPmR', myTime, myThid )
231         IF ( balanceQnet )        IF ( balanceQnet )
232       &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,       &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,
233       &        'Qnet ', myTime, myThid )       &        'Qnet ', myTime, myThid )
234  #endif /* ALLOW_BALANCE_FLUXES */  #endif /* ALLOW_BALANCE_FLUXES */
# Line 196  CHPF$ INDEPENDENT Line 255  CHPF$ INDEPENDENT
255            itdkey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
256       &                      + act3*max1*max2       &                      + act3*max1*max2
257       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
258    #else  /* ALLOW_AUTODIFF_TAMC */
259    C     if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
260    C     and all vertical mixing schemes, but keep OBCS_CALC
261            IF ( fluidIsWater ) THEN
262  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
263    
264  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 208  C     uninitialised but inert locations. Line 271  C     uninitialised but inert locations.
271           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
272            rhok   (i,j)   = 0. _d 0            rhok   (i,j)   = 0. _d 0
273            rhoKM1 (i,j)   = 0. _d 0            rhoKM1 (i,j)   = 0. _d 0
274              rhoKP1 (i,j)   = 0. _d 0
275           ENDDO           ENDDO
276          ENDDO          ENDDO
277    
# Line 264  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_ Line 328  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_
328  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
329    
330  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
331          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
332       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
333  #endif  #endif
334    
# Line 281  cph Needed for rhok, rhokm1, in the case Line 345  cph Needed for rhok, rhokm1, in the case
345    
346  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
347  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
348            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
349       &                   .OR. doDiagsRho.GE.1 ) THEN       &                   .OR. doDiagsRho.GE.1 ) THEN
350  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
351              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
352       &       CALL DEBUG_CALL('FIND_RHO',myThid)       &       CALL DEBUG_CALL('FIND_RHO',myThid)
353  #endif  #endif
354  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 310  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 373  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
373       I        myThid )       I        myThid )
374              ENDIF              ENDIF
375  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
376              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
377       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
378  #endif  #endif
379    cph Avoid variable aliasing for adjoint !!!
380                DO j=jMin,jMax
381                 DO i=iMin,iMax
382                  rhoKP1(i,j) = rhoK(i,j)
383                 ENDDO
384                ENDDO
385              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
386       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
387       I             rhoK, rhoKm1, rhoK,       I             rhoK, rhoKm1, rhoKp1,
388       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
389       I             myThid )       I             myThid )
390            ENDIF            ENDIF
# Line 330  C--       Implicit Vertical Diffusion fo Line 399  C--       Implicit Vertical Diffusion fo
399  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
400            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
401  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
402              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
403       &       CALL DEBUG_CALL('CALC_IVDC',myThid)       &       CALL DEBUG_CALL('CALC_IVDC',myThid)
404  #endif  #endif
405              CALL CALC_IVDC(              CALL CALC_IVDC(
# Line 351  C--     end of diagnostic k loop (Nr:1) Line 420  C--     end of diagnostic k loop (Nr:1)
420          ENDDO          ENDDO
421    
422  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
 c       IF ( useDiagnostics .AND.  
 c    &       (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN  
423          IF ( doDiagsRho.GE.1 ) THEN          IF ( doDiagsRho.GE.1 ) THEN
424            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
425       &         2, bi, bj, myThid)       &         2, bi, bj, myThid)
426          ENDIF          ENDIF
427  #endif  #endif
428    
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
 #ifdef ALLOW_DEBUG  
           IF ( debugLevel .GE. debLevB )  
      &     CALL DEBUG_CALL('OBCS_CALC',myThid)  
 #endif  
           CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 #ifndef ALLOW_AUTODIFF_TAMC  
         IF ( fluidIsWater ) THEN  
 #endif  
429  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
430  C       relaxation terms, etc.  C       relaxation terms, etc.
431  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
432          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
433       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
434  #endif  #endif
435  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 395  CADJ STORE recip_hFacC(:,:,:,bi,bj) Line 446  CADJ STORE recip_hFacC(:,:,:,bi,bj)
446  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
447  # endif  # endif
448  #endif  #endif
449           CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
450       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
451       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
 #ifndef ALLOW_AUTODIFF_TAMC  
         ENDIF  
 #endif  
452  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
453  # ifdef EXACT_CONSERV  # ifdef EXACT_CONSERV
454  cph-test  cph-test
# Line 440  CADJ STORE sigmaR(:,:,:)        = comlev Line 488  CADJ STORE sigmaR(:,:,:)        = comlev
488  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
489          IF (useGMRedi) THEN          IF (useGMRedi) THEN
490  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
491            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
492       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
493  #endif  #endif
494            CALL GMREDI_CALC_TENSOR(            CALL GMREDI_CALC_TENSOR(
# Line 462  C--     Calculate iso-neutral slopes for Line 510  C--     Calculate iso-neutral slopes for
510  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
511          IF (useKPP) THEN          IF (useKPP) THEN
512  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
513            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
514       &     CALL DEBUG_CALL('KPP_CALC',myThid)       &     CALL DEBUG_CALL('KPP_CALC',myThid)
515  #endif  #endif
516            CALL KPP_CALC(            CALL KPP_CALC(
# Line 480  C--     Compute KPP mixing coefficients Line 528  C--     Compute KPP mixing coefficients
528  C--     Compute PP81 mixing coefficients  C--     Compute PP81 mixing coefficients
529          IF (usePP81) THEN          IF (usePP81) THEN
530  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
531            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
532       &     CALL DEBUG_CALL('PP81_CALC',myThid)       &     CALL DEBUG_CALL('PP81_CALC',myThid)
533  #endif  #endif
534            CALL PP81_CALC(            CALL PP81_CALC(
# Line 492  C--     Compute PP81 mixing coefficients Line 540  C--     Compute PP81 mixing coefficients
540  C--     Compute MY82 mixing coefficients  C--     Compute MY82 mixing coefficients
541          IF (useMY82) THEN          IF (useMY82) THEN
542  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
543            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
544       &     CALL DEBUG_CALL('MY82_CALC',myThid)       &     CALL DEBUG_CALL('MY82_CALC',myThid)
545  #endif  #endif
546            CALL MY82_CALC(            CALL MY82_CALC(
# Line 504  C--     Compute MY82 mixing coefficients Line 552  C--     Compute MY82 mixing coefficients
552  C--     Compute GGL90 mixing coefficients  C--     Compute GGL90 mixing coefficients
553          IF (useGGL90) THEN          IF (useGGL90) THEN
554  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
555            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
556       &     CALL DEBUG_CALL('GGL90_CALC',myThid)       &     CALL DEBUG_CALL('GGL90_CALC',myThid)
557  #endif  #endif
558            CALL GGL90_CALC(            CALL GGL90_CALC(
# Line 513  C--     Compute GGL90 mixing coefficient Line 561  C--     Compute GGL90 mixing coefficient
561  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
562    
563  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
564          IF ( taveFreq.GT. 0. _d 0 .AND. fluidIsWater ) THEN          IF ( taveFreq.GT. 0. _d 0 ) THEN
565            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
566          ENDIF          ENDIF
567          IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN          IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
# Line 522  C--     Compute GGL90 mixing coefficient Line 570  C--     Compute GGL90 mixing coefficient
570          ENDIF          ENDIF
571  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
572    
573    #ifndef ALLOW_AUTODIFF_TAMC
574    C---  if fluid Is Water: end
575            ENDIF
576    #endif
577    
578    #ifdef  ALLOW_OBCS
579    C--     Calculate future values on open boundaries
580            IF (useOBCS) THEN
581    #ifdef ALLOW_DEBUG
582              IF ( debugLevel .GE. debLevB )
583         &     CALL DEBUG_CALL('OBCS_CALC',myThid)
584    #endif
585              CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
586         I            uVel, vVel, wVel, theta, salt,
587         I            myThid )
588            ENDIF
589    #endif  /* ALLOW_OBCS */
590    
591  C--   end bi,bj loops.  C--   end bi,bj loops.
592         ENDDO         ENDDO
593        ENDDO        ENDDO
# Line 537  C--   end bi,bj loops. Line 603  C--   end bi,bj loops.
603  #endif  #endif
604    
605  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
606           IF ( debugLevel .GE. debLevB )        IF ( debugLevel .GE. debLevB )
607       &    CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)       &     CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
608  #endif  #endif
609    
610        RETURN        RETURN
Legend:
Removed from v.1.28  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22