/[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.30 by jscott, Wed Sep 6 15:30:25 2006 UTC revision 1.46 by heimbach, Fri Jun 1 22:20:02 2007 UTC
# Line 64  C     == Global variables === Line 64  C     == Global variables ===
64  # endif  # endif
65  # ifdef ALLOW_EXF  # ifdef ALLOW_EXF
66  #  include "ctrl.h"  #  include "ctrl.h"
67  #  include "exf_fields.h"  #  include "EXF_FIELDS.h"
 #  include "exf_clim_fields.h"  
68  #  ifdef ALLOW_BULKFORMULAE  #  ifdef ALLOW_BULKFORMULAE
69  #   include "exf_constants.h"  #   include "EXF_CONSTANTS.h"
70  #  endif  #  endif
71  # endif  # endif
72  # ifdef ALLOW_SEAICE  # ifdef ALLOW_SEAICE
# Line 91  C     iMin, iMax    :: Ranges and sub-bl Line 90  C     iMin, iMax    :: Ranges and sub-bl
90  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
91  C     bi, bj        :: tile indices  C     bi, bj        :: tile indices
92  C     i,j,k         :: loop indices  C     i,j,k         :: loop indices
93          _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 117  C--   dummy statement to end declaration Line 117  C--   dummy statement to end declaration
117        IF ( debugLevel .GE. debLevB )        IF ( debugLevel .GE. debLevB )
118       &     CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)       &     CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
119  #endif  #endif
120    
121        doDiagsRho = 0        doDiagsRho = 0
122  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
123        IF ( useDiagnostics .AND. fluidIsWater ) THEN        IF ( useDiagnostics .AND. fluidIsWater ) THEN
# Line 143  C     under sea-ice may not be "exactly" Line 143  C     under sea-ice may not be "exactly"
143  C     theta is dumped or time-averaged.  C     theta is dumped or time-averaged.
144        IF ( useSEAICE ) THEN        IF ( useSEAICE ) THEN
145  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
146  CADJ STORE aqh,precip,swdown   = comlev1, key = ikey_dynamics  CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics
147    CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics
148  CADJ STORE theta               = comlev1, key = ikey_dynamics  CADJ STORE theta               = comlev1, key = ikey_dynamics
149  # ifdef SEAICE_ALLOW_DYNAMICS  cph# ifdef EXF_READ_EVAP
150    CADJ STORE evap                = comlev1, key = ikey_dynamics
151    cph# endif
152    cph# ifdef SEAICE_ALLOW_DYNAMICS
153  CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics  CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics
154    CADJ STORE uice,vice           = comlev1, key = ikey_dynamics
155    CADJ STORE zeta, eta           = comlev1, key = ikey_dynamics
156    cph# endif
157    # ifdef SEAICE_CGRID
158    CADJ STORE fu, fv              = comlev1, key = ikey_dynamics
159    CADJ STORE seaicemasku         = comlev1, key = ikey_dynamics
160    CADJ STORE seaicemaskv         = comlev1, key = ikey_dynamics
161    #  ifdef ATMOSPHERIC_LOADING
162    CADJ STORE siceload            = comlev1, key = ikey_dynamics
163    #  endif
164    # endif
165    # ifdef SEAICE_ALLOW_EVP
166    CADJ STORE dwatn           = comlev1, key = ikey_dynamics
167    CADJ STORE seaice_sigma1   = comlev1, key = ikey_dynamics
168    CADJ STORE seaice_sigma2   = comlev1, key = ikey_dynamics
169    CADJ STORE seaice_sigma12  = comlev1, key = ikey_dynamics
170  # endif  # endif
171  #endif  #endif
172  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
# Line 205  CADJ STORE theta = comlev1, key = ikey_d Line 225  CADJ STORE theta = comlev1, key = ikey_d
225  #ifdef ALLOW_OCN_COMPON_INTERF  #ifdef ALLOW_OCN_COMPON_INTERF
226  C--    Apply imported data (from coupled interface) to forcing fields  C--    Apply imported data (from coupled interface) to forcing fields
227  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 ?)
228         IF ( useCoupler ) THEN        IF ( useCoupler ) THEN
229           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
230         ENDIF        ENDIF
231  #endif /* ALLOW_OCN_COMPON_INTERF */  #endif /* ALLOW_OCN_COMPON_INTERF */
232    
233  #ifdef ALLOW_BALANCE_FLUXES  #ifdef ALLOW_BALANCE_FLUXES
234  C     balance fluxes  C     balance fluxes
235         IF ( balanceEmPmR )        IF ( balanceEmPmR )
236       &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,       &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,
237       &        'EmPmR', myTime, myThid )       &        'EmPmR', myTime, myThid )
238         IF ( balanceQnet )        IF ( balanceQnet )
239       &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,       &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,
240       &        'Qnet ', myTime, myThid )       &        'Qnet ', myTime, myThid )
241  #endif /* ALLOW_BALANCE_FLUXES */  #endif /* ALLOW_BALANCE_FLUXES */
# Line 242  CHPF$ INDEPENDENT Line 262  CHPF$ INDEPENDENT
262            itdkey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
263       &                      + act3*max1*max2       &                      + act3*max1*max2
264       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
265    #else  /* ALLOW_AUTODIFF_TAMC */
266    C     if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
267    C     and all vertical mixing schemes, but keep OBCS_CALC
268            IF ( fluidIsWater ) THEN
269  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
270    
271  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 254  C     uninitialised but inert locations. Line 278  C     uninitialised but inert locations.
278           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
279            rhok   (i,j)   = 0. _d 0            rhok   (i,j)   = 0. _d 0
280            rhoKM1 (i,j)   = 0. _d 0            rhoKM1 (i,j)   = 0. _d 0
281              rhoKP1 (i,j)   = 0. _d 0
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284    
# Line 288  cph although some of these are re-initia Line 313  cph although some of these are re-initia
313             VisbeckK(i,j,bi,bj)   = 0. _d 0             VisbeckK(i,j,bi,bj)   = 0. _d 0
314  #  endif  #  endif
315  # endif /* ALLOW_GMREDI */  # endif /* ALLOW_GMREDI */
316    # ifdef ALLOW_KPP
317               KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0
318               KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0
319    # endif /* ALLOW_KPP */
320  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
321            ENDDO            ENDDO
322           ENDDO           ENDDO
# Line 310  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_ Line 339  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_
339  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
340    
341  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
342          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
343       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
344  #endif  #endif
345    
# Line 327  cph Needed for rhok, rhokm1, in the case Line 356  cph Needed for rhok, rhokm1, in the case
356    
357  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
358  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  
359            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
360       &                   .OR. doDiagsRho.GE.1 ) THEN       &                   .OR. doDiagsRho.GE.1 ) THEN
361  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
362              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
363       &       CALL DEBUG_CALL('FIND_RHO',myThid)       &       CALL DEBUG_CALL('FIND_RHO',myThid)
364  #endif  #endif
365  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 356  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 384  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
384       I        myThid )       I        myThid )
385              ENDIF              ENDIF
386  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
387              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
388       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
389  #endif  #endif
390    cph Avoid variable aliasing for adjoint !!!
391                DO j=jMin,jMax
392                 DO i=iMin,iMax
393                  rhoKP1(i,j) = rhoK(i,j)
394                 ENDDO
395                ENDDO
396              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
397       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
398       I             rhoK, rhoKm1, rhoK,       I             rhoK, rhoKm1, rhoKp1,
399       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
400       I             myThid )       I             myThid )
401            ENDIF            ENDIF
# Line 376  C--       Implicit Vertical Diffusion fo Line 410  C--       Implicit Vertical Diffusion fo
410  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
411            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
412  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
413              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
414       &       CALL DEBUG_CALL('CALC_IVDC',myThid)       &       CALL DEBUG_CALL('CALC_IVDC',myThid)
415  #endif  #endif
416              CALL CALC_IVDC(              CALL CALC_IVDC(
# Line 397  C--     end of diagnostic k loop (Nr:1) Line 431  C--     end of diagnostic k loop (Nr:1)
431          ENDDO          ENDDO
432    
433  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
 c       IF ( useDiagnostics .AND.  
 c    &       (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN  
434          IF ( doDiagsRho.GE.1 ) THEN          IF ( doDiagsRho.GE.1 ) THEN
435            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
436       &         2, bi, bj, myThid)       &         2, bi, bj, myThid)
437          ENDIF          ENDIF
438  #endif  #endif
439    
 #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  
440  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
441  C       relaxation terms, etc.  C       relaxation terms, etc.
442  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
443          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
444       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
445  #endif  #endif
446  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 441  CADJ STORE recip_hFacC(:,:,:,bi,bj) Line 457  CADJ STORE recip_hFacC(:,:,:,bi,bj)
457  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
458  # endif  # endif
459  #endif  #endif
460           CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
461       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
462       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
 #ifndef ALLOW_AUTODIFF_TAMC  
         ENDIF  
 #endif  
463  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
464  # ifdef EXACT_CONSERV  # ifdef EXACT_CONSERV
465  cph-test  cph-test
# Line 486  CADJ STORE sigmaR(:,:,:)        = comlev Line 499  CADJ STORE sigmaR(:,:,:)        = comlev
499  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
500          IF (useGMRedi) THEN          IF (useGMRedi) THEN
501  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
502            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
503       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
504  #endif  #endif
505            CALL GMREDI_CALC_TENSOR(            CALL GMREDI_CALC_TENSOR(
# Line 508  C--     Calculate iso-neutral slopes for Line 521  C--     Calculate iso-neutral slopes for
521  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
522          IF (useKPP) THEN          IF (useKPP) THEN
523  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
524            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
525       &     CALL DEBUG_CALL('KPP_CALC',myThid)       &     CALL DEBUG_CALL('KPP_CALC',myThid)
526  #endif  #endif
527            CALL KPP_CALC(            CALL KPP_CALC(
528       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
529  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
530          ELSE          ELSE
531            CALL KPP_CALC_DUMMY(            CALL KPP_CALC_DUMMY(
532       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
533  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
534          ENDIF          ENDIF
535    
# Line 526  C--     Compute KPP mixing coefficients Line 539  C--     Compute KPP mixing coefficients
539  C--     Compute PP81 mixing coefficients  C--     Compute PP81 mixing coefficients
540          IF (usePP81) THEN          IF (usePP81) THEN
541  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
542            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
543       &     CALL DEBUG_CALL('PP81_CALC',myThid)       &     CALL DEBUG_CALL('PP81_CALC',myThid)
544  #endif  #endif
545            CALL PP81_CALC(            CALL PP81_CALC(
# Line 538  C--     Compute PP81 mixing coefficients Line 551  C--     Compute PP81 mixing coefficients
551  C--     Compute MY82 mixing coefficients  C--     Compute MY82 mixing coefficients
552          IF (useMY82) THEN          IF (useMY82) THEN
553  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
554            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
555       &     CALL DEBUG_CALL('MY82_CALC',myThid)       &     CALL DEBUG_CALL('MY82_CALC',myThid)
556  #endif  #endif
557            CALL MY82_CALC(            CALL MY82_CALC(
# Line 550  C--     Compute MY82 mixing coefficients Line 563  C--     Compute MY82 mixing coefficients
563  C--     Compute GGL90 mixing coefficients  C--     Compute GGL90 mixing coefficients
564          IF (useGGL90) THEN          IF (useGGL90) THEN
565  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
566            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
567       &     CALL DEBUG_CALL('GGL90_CALC',myThid)       &     CALL DEBUG_CALL('GGL90_CALC',myThid)
568  #endif  #endif
569            CALL GGL90_CALC(            CALL GGL90_CALC(
# Line 559  C--     Compute GGL90 mixing coefficient Line 572  C--     Compute GGL90 mixing coefficient
572  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
573    
574  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
575          IF ( taveFreq.GT. 0. _d 0 .AND. fluidIsWater ) THEN          IF ( taveFreq.GT. 0. _d 0 ) THEN
576            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
577          ENDIF          ENDIF
578          IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN          IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
# Line 568  C--     Compute GGL90 mixing coefficient Line 581  C--     Compute GGL90 mixing coefficient
581          ENDIF          ENDIF
582  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
583    
584    #ifndef ALLOW_AUTODIFF_TAMC
585    C---  if fluid Is Water: end
586            ENDIF
587    #endif
588    
589    #ifdef  ALLOW_OBCS
590    C--     Calculate future values on open boundaries
591            IF (useOBCS) THEN
592    #ifdef ALLOW_DEBUG
593              IF ( debugLevel .GE. debLevB )
594         &     CALL DEBUG_CALL('OBCS_CALC',myThid)
595    #endif
596              CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
597         I            uVel, vVel, wVel, theta, salt,
598         I            myThid )
599            ENDIF
600    #endif  /* ALLOW_OBCS */
601    
602  C--   end bi,bj loops.  C--   end bi,bj loops.
603         ENDDO         ENDDO
604        ENDDO        ENDDO
605    
606    #ifdef  ALLOW_KPP
607          IF (useKPP) THEN
608            CALL KPP_DO_EXCH( myThid )
609          ENDIF
610    #endif  /* ALLOW_KPP */
611    
612  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
613        IF ( fluidIsWater .AND. useDiagnostics ) THEN        IF ( fluidIsWater .AND. useDiagnostics ) THEN
614          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
Legend:
Removed from v.1.30  
changed lines
  Added in v.1.46
  ViewVC Help
Powered by ViewVC 1.1.22