/[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.66 by gforget, Fri May 30 21:10:31 2008 UTC revision 1.74 by jmc, Fri Oct 3 02:26:49 2008 UTC
# Line 36  C     == Global variables === Line 36  C     == Global variables ===
36  #include "SIZE.h"  #include "SIZE.h"
37  #include "EEPARAMS.h"  #include "EEPARAMS.h"
38  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
39  #include "GRID.h"  #include "GRID.h"
40    #include "DYNVARS.h"
41  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
42  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
43  #endif  #endif
# Line 90  C     bi, bj        :: tile indices Line 90  C     bi, bj        :: tile indices
90  C     i,j,k         :: loop indices  C     i,j,k         :: loop indices
91        _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL rhoKm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhoKm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhoK    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
93        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 119  C--   dummy statement to end declaration Line 118  C--   dummy statement to end declaration
118        doDiagsRho = 0        doDiagsRho = 0
119  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
120        IF ( useDiagnostics .AND. fluidIsWater ) THEN        IF ( useDiagnostics .AND. fluidIsWater ) THEN
121          IF ( DIAGNOSTICS_IS_ON('RHOANOSQ',myThid) .OR.          IF ( DIAGNOSTICS_IS_ON('WRHOMASS',myThid) )
122       &       DIAGNOSTICS_IS_ON('URHOMASS',myThid) .OR.       &       doDiagsRho = doDiagsRho + 1
123       &       DIAGNOSTICS_IS_ON('VRHOMASS',myThid) .OR.          IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) )
124       &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) .OR.       &       doDiagsRho = doDiagsRho + 2
125       &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) ) doDiagsRho = 2          IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
126          IF ( doDiagsRho.EQ.0 .AND.       &       doDiagsRho = doDiagsRho + 4
      &       DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) ) doDiagsRho = 1  
         IF ( doDiagsRho.EQ.0 .AND.  
      &       DIAGNOSTICS_IS_ON('DRHODR  ',myThid) ) doDiagsRho = 1  
127        ENDIF        ENDIF
128  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
129    
130    
131  #ifdef ALLOW_SEAICE  #ifdef ALLOW_SEAICE
132        IF ( useSEAICE ) THEN        IF ( useSEAICE ) THEN
133  # ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
# Line 218  C--   Freeze water at the surface Line 215  C--   Freeze water at the surface
215  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
216  CADJ STORE theta = comlev1, key = ikey_dynamics  CADJ STORE theta = comlev1, key = ikey_dynamics
217  #endif  #endif
218        IF ( allowFreezing        IF ( allowFreezing ) THEN
      &                   .AND. .NOT. useSEAICE  
      &                   .AND. .NOT. useThSIce ) THEN  
219          CALL FREEZE_SURFACE(  myTime, myIter, myThid )          CALL FREEZE_SURFACE(  myTime, myIter, myThid )
220        ENDIF        ENDIF
221    
# Line 276  C     just ensure that all memory refere Line 271  C     just ensure that all memory refere
271  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
272  C     uninitialised but inert locations.  C     uninitialised but inert locations.
273    
274    #ifdef ALLOW_AUTODIFF_TAMC
275          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
276           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
           rhoK   (i,j)   = 0. _d 0  
277            rhoKm1 (i,j)   = 0. _d 0            rhoKm1 (i,j)   = 0. _d 0
278            rhoKp1 (i,j)   = 0. _d 0            rhoKp1 (i,j)   = 0. _d 0
279           ENDDO           ENDDO
280          ENDDO          ENDDO
281    #endif /* ALLOW_AUTODIFF_TAMC */
282    
283          DO k=1,Nr          DO k=1,Nr
284           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
285            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
286  C This is currently also used by IVDC and Diagnostics  C This is currently used by GMRedi, IVDC, MXL-depth  and Diagnostics
287             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
288             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
289             sigmaR(i,j,k) = 0. _d 0             sigmaR(i,j,k) = 0. _d 0
290  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
291  cph all the following init. are necessary for TAF  cph all the following init. are necessary for TAF
292  cph although some of these are re-initialised later.  cph although some of these are re-initialised later.
293    c          rhoInSitu(i,j,k,bi,bj) = 0.
294             IVDConvCount(i,j,k,bi,bj) = 0.             IVDConvCount(i,j,k,bi,bj) = 0.
295  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
296             Kwx(i,j,k,bi,bj)  = 0. _d 0             Kwx(i,j,k,bi,bj)  = 0. _d 0
# Line 351  C--     Start of diagnostic loop Line 348  C--     Start of diagnostic loop
348  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
349  C? Patrick, is this formula correct now that we change the loop range?  C? Patrick, is this formula correct now that we change the loop range?
350  C? Do we still need this?  C? Do we still need this?
351  cph kkey formula corrected.  cph kkey formula corrected.
352  cph Needed for rhoK, rhoKm1, in the case useGMREDI.  cph Needed for rhoK, rhoKm1, in the case useGMREDI.
353           kkey = (itdkey-1)*Nr + k            kkey = (itdkey-1)*Nr + k
354  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
355    
356  C--       Calculate gradients of potential density for isoneutral  C--   Always compute density (stored in common block) here; even when it is not
357  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C     needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)  
      &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN  
358  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
359              IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
360       &       CALL DEBUG_CALL('FIND_RHO',myThid)       &       CALL DEBUG_CALL('FIND_RHO_2D',myThid)
361  #endif  #endif
362  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
363              IF ( fluidIsWater ) THEN
364  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
365  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
366    #endif /* ALLOW_AUTODIFF_TAMC */
367    #ifdef ALLOW_DOWN_SLOPE
368              IF ( useDOWN_SLOPE ) THEN
369                CALL DWNSLP_CALC_RHO(
370         I                  theta, salt,
371         O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
372         I                  k, bi, bj, myTime, myIter, myThid )
373              ELSE
374    #endif /* ALLOW_DOWN_SLOPE */
375                CALL FIND_RHO_2D(
376         I                iMin, iMax, jMin, jMax, k,
377         I                theta(1-OLx,1-OLy,k,bi,bj),
378         I                salt (1-OLx,1-OLy,k,bi,bj),
379         O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
380         I                k, bi, bj, myThid )
381    #ifdef ALLOW_DOWN_SLOPE
382              ENDIF
383    #endif /* ALLOW_DOWN_SLOPE */
384    #ifdef ALLOW_AUTODIFF_TAMC
385              ELSE
386    C-        fluid is not water:
387               DO j=1-OLy,sNy+OLy
388                DO i=1-OLx,sNx+OLx
389                  rhoInSitu(i,j,k,bi,bj) = 0.
390                ENDDO
391               ENDDO
392              ENDIF
393  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
394    
395    C--       Calculate gradients of potential density for isoneutral
396    C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
397              IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
398         &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
399              IF (k.GT.1) THEN              IF (k.GT.1) THEN
400  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
401  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
402  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
403    CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey, byte=isbyte
404  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
405               CALL FIND_RHO(               CALL FIND_RHO_2D(
406       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,       I                 iMin, iMax, jMin, jMax, k,
407       I        theta, salt,       I                 theta(1-OLx,1-OLy,k-1,bi,bj),
408       O        rhoKm1,       I                 salt (1-OLx,1-OLy,k-1,bi,bj),
409       I        myThid )       O                 rhoKm1,
410         I                 k-1, bi, bj, myThid )
411              ENDIF              ENDIF
412  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
413              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
# Line 392  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 416  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
416  cph Avoid variable aliasing for adjoint !!!  cph Avoid variable aliasing for adjoint !!!
417              DO j=jMin,jMax              DO j=jMin,jMax
418               DO i=iMin,iMax               DO i=iMin,iMax
419                rhoKp1(i,j) = rhoK(i,j)                rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
420               ENDDO               ENDDO
421              ENDDO              ENDDO
422              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
423       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
424       I             rhoK, rhoKm1, rhoKp1,       I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
425       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
426       I             myThid )       I             myThid )
427  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
428  #ifdef GMREDI_WITH_STABLE_ADJOINT  #ifdef GMREDI_WITH_STABLE_ADJOINT
429  cgf zero out adjoint fields to stabilize pkg/gmredi adjoint  cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
430  cgf -> cuts adjoint dependency from slope to state  cgf -> cuts adjoint dependency from slope to state
431        CALL ZERO_ADJ_3D( bi, bj, sigmaX, myThid)              CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
432        CALL ZERO_ADJ_3D( bi, bj, sigmaY, myThid)              CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
433        CALL ZERO_ADJ_3D( bi, bj, sigmaR, myThid)              CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
434  #endif  #endif
435  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
436            ENDIF            ENDIF
# Line 420  c ==> should use sigmaR !!! Line 444  c ==> should use sigmaR !!!
444  #endif  #endif
445              CALL CALC_IVDC(              CALL CALC_IVDC(
446       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
447       I        rhoKm1, rhoK,       I        rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
448       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
449            ENDIF            ENDIF
450    
451  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
452            IF ( doDiagsRho.GE.2 ) THEN            IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
453              CALL DIAGS_RHO( k, bi, bj,              CALL DIAGS_RHO_L( k, bi, bj,
454       I                      rhoK, rhoKm1,       I                        rhoInSitu(1-OLx,1-OLy,k,bi,bj),
455       I                      myTime, myIter, myThid)       I                        rhoKm1, wVel,
456         I                        myTime, myIter, myThid )
457            ENDIF            ENDIF
458  #endif  #endif
459    
# Line 436  C--     end of diagnostic k loop (Nr:1) Line 461  C--     end of diagnostic k loop (Nr:1)
461          ENDDO          ENDDO
462    
463  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
464  CADJ STORE IVDConvCount(:,:,:,bi,bj)  CADJ STORE IVDConvCount(:,:,:,bi,bj)
465  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
466  #endif  #endif
467    
468  C--     Diagnose Mixed Layer Depth:  C--     Diagnose Mixed Layer Depth:
469          IF ( useGMRedi .OR. doDiagsRho.GE.1 ) THEN          IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
470            CALL CALC_OCE_MXLAYER( rhoK, sigmaR,            CALL CALC_OCE_MXLAYER(
471       &              bi, bj, myTime, myIter, myThid )       I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
472         I              bi, bj, myTime, myIter, myThid )
473          ENDIF          ENDIF
474    
475  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
476          IF ( useSALT_PLUME ) THEN          IF ( useSALT_PLUME ) THEN
477            CALL SALT_PLUME_CALC_DEPTH( rhoK, sigmaR,            CALL SALT_PLUME_CALC_DEPTH(
478       &              bi, bj, myTime, myIter, myThid )       I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
479         I              bi, bj, myTime, myIter, myThid )
480          ENDIF          ENDIF
481  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
482    
483  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
484          IF ( doDiagsRho.GE.1 ) THEN          IF ( MOD(doDiagsRho,4).GE.2 ) THEN
485            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
486       &         2, bi, bj, myThid)       &         2, bi, bj, myThid)
487          ENDIF          ENDIF
# Line 569  C--     Compute GGL90 mixing coefficient Line 596  C--     Compute GGL90 mixing coefficient
596          ENDIF          ENDIF
597  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
598    
599  #ifdef  ALLOW_GMREDI  #ifdef ALLOW_GMREDI
600  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
601  # ifndef GM_EXCLUDE_CLIPPING  # ifndef GM_EXCLUDE_CLIPPING
602  cph storing here is needed only for one GMREDI_OPTIONS:  cph storing here is needed only for one GMREDI_OPTIONS:
# Line 589  C--     Calculate iso-neutral slopes for Line 616  C--     Calculate iso-neutral slopes for
616       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
617  #endif  #endif
618            CALL GMREDI_CALC_TENSOR(            CALL GMREDI_CALC_TENSOR(
 c    I             bi, bj, iMin, iMax, jMin, jMax,  
 c    I             sigmaX, sigmaY, sigmaR,  
 c    I             myThid )  
619       I             iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
620       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
621       I             bi, bj, myTime, myIter, myThid )       I             bi, bj, myTime, myIter, myThid )
622  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
623          ELSE          ELSE
624            CALL GMREDI_CALC_TENSOR_DUMMY(            CALL GMREDI_CALC_TENSOR_DUMMY(
 c    I             bi, bj, iMin, iMax, jMin, jMax,  
 c    I             sigmaX, sigmaY, sigmaR,  
 c    I             myThid )  
625       I             iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
626       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
627       I             bi, bj, myTime, myIter, myThid )       I             bi, bj, myTime, myIter, myThid )
628  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
629          ENDIF          ENDIF
630  #endif  /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
631    
632    #ifdef ALLOW_DOWN_SLOPE
633            IF ( useDOWN_SLOPE ) THEN
634    C--     Calculate Downsloping Flow for Down_Slope parameterization
635             IF ( usingPCoords ) THEN
636              CALL DWNSLP_CALC_FLOW(
637         I                bi, bj, kSurfC, rhoInSitu,
638         I                myTime, myIter, myThid )
639             ELSE
640              CALL DWNSLP_CALC_FLOW(
641         I                bi, bj, kLowC, rhoInSitu,
642         I                myTime, myIter, myThid )
643             ENDIF
644            ENDIF
645    #endif /* ALLOW_DOWN_SLOPE */
646    
647  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
648  C---  if fluid Is Water: end  C---  if fluid Is Water: end
# Line 638  C--   end bi,bj loops. Line 674  C--   end bi,bj loops.
674    
675  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
676        IF ( fluidIsWater .AND. useDiagnostics ) THEN        IF ( fluidIsWater .AND. useDiagnostics ) THEN
677            CALL DIAGS_RHO_G(
678         I                    rhoInSitu, uVel, vVel,
679         I                    myTime, myIter, myThid )
680          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
681        ENDIF        ENDIF
682        IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN        IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
683          CALL DIAGNOSTICS_FILL( IVDConvCount,'CONVADJ ',          CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
684       &                         0, Nr, 0, 1, 1, myThid )       &                               0, Nr, 0, 1, 1, myThid )
685        ENDIF        ENDIF
686  #endif  #endif
687    
Legend:
Removed from v.1.66  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22