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

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

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

revision 1.41 by jmc, Wed Apr 11 04:02:05 2012 UTC revision 1.42 by jmc, Tue Dec 18 01:16:53 2012 UTC
# Line 78  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy Line 78  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy
78  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
79  C     == Local variables ==  C     == Local variables ==
80        INTEGER i,j        INTEGER i,j
       _RL zero, one, half  
81        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL dRlocM,dRlocP, ddRloc, locAlpha        _RL dRlocM,dRlocP, ddRloc, locAlpha
83        _RL ddPIm, ddPIp, rec_dRm, rec_dRp        _RL ddPIm, ddPIp, rec_dRm, rec_dRp
84        _RL surfPhiFac        _RL surfPhiFac
       PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )  
85        LOGICAL useDiagPhiRlow, addSurfPhiAnom        LOGICAL useDiagPhiRlow, addSurfPhiAnom
86  CEOP  CEOP
87        useDiagPhiRlow = .TRUE.        useDiagPhiRlow = .TRUE.
# Line 184  C    but might be slower (--> keep origi Line 182  C    but might be slower (--> keep origi
182  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
183    
184  #ifdef ALLOW_MOM_COMMON  #ifdef ALLOW_MOM_COMMON
185  C Quasi-hydrostatic terms are added in as if they modify the buoyancy  C--     Quasi-hydrostatic terms are added in as if they modify the buoyancy
186          IF (quasiHydrostatic) THEN          IF (quasiHydrostatic) THEN
187           CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)           CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
188          ENDIF          ENDIF
# Line 215  C           conserve KE+PE exactly even Line 213  C           conserve KE+PE exactly even
213           DO j=jMin,jMax           DO j=jMin,jMax
214            DO i=iMin,iMax            DO i=iMin,iMax
215              phiHydC(i,j) = phiHydF(i,j)              phiHydC(i,j) = phiHydF(i,j)
216       &            + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst       &          + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
217              phiHydF(i,j) = phiHydF(i,j)              phiHydF(i,j) = phiHydF(i,j)
218       &                 + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst       &                 + drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
219            ENDDO            ENDDO
# Line 231  C           conserve KE+PE exactly even Line 229  C           conserve KE+PE exactly even
229              phiHydC(i,j) = ddRloc*gravity*alphaRho(i,j)*recip_rhoConst              phiHydC(i,j) = ddRloc*gravity*alphaRho(i,j)*recip_rhoConst
230             ELSE             ELSE
231              phiHydC(i,j) = phiHydF(i,j)              phiHydC(i,j) = phiHydF(i,j)
232       &              + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst       &            + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
233             ENDIF             ENDIF
234              phiHydF(i,j) = phiHydC(i,j)              phiHydF(i,j) = phiHydC(i,j)
235       &              + half*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst       &            + halfRL*drF(k)*gravity*alphaRho(i,j)*recip_rhoConst
236            ENDDO            ENDDO
237           ENDDO           ENDDO
238          ENDIF          ENDIF
# Line 245  C  --  Finite Difference Form Line 243  C  --  Finite Difference Form
243  C---------- This discretization is the "energy conserving" form  C---------- This discretization is the "energy conserving" form
244  C           which has been used since at least Adcroft et al., MWR 1997  C           which has been used since at least Adcroft et al., MWR 1997
245    
246          dRlocM=half*drC(k)          dRlocM = halfRL*drC(k)
247          IF (k.EQ.1) dRlocM=rF(k)-rC(k)          IF (k.EQ.1) dRlocM=rF(k)-rC(k)
248          IF (k.EQ.Nr) THEN          IF (k.EQ.Nr) THEN
249            dRlocP=rC(k)-rF(k+1)            dRlocP=rC(k)-rF(k+1)
250          ELSE          ELSE
251            dRlocP=half*drC(k+1)            dRlocP=halfRL*drC(k+1)
252          ENDIF          ENDIF
253          IF ( uniformFreeSurfLev ) THEN          IF ( uniformFreeSurfLev ) THEN
254           DO j=jMin,jMax           DO j=jMin,jMax
# Line 262  C           which has been used since at Line 260  C           which has been used since at
260            ENDDO            ENDDO
261           ENDDO           ENDDO
262          ELSE          ELSE
263           rec_dRm = one/(rF(k)-rC(k))           rec_dRm = oneRL/(rF(k)-rC(k))
264           rec_dRp = one/(rC(k)-rF(k+1))           rec_dRp = oneRL/(rC(k)-rF(k+1))
265           DO j=jMin,jMax           DO j=jMin,jMax
266            DO i=iMin,iMax            DO i=iMin,iMax
267             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
# Line 271  C           which has been used since at Line 269  C           which has been used since at
269  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
270              ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)              ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
271  #endif  #endif
272              phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM              phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
273       &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP       &                     +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
274       &                     )*gravity*alphaRho(i,j)*recip_rhoConst       &                    )*gravity*alphaRho(i,j)*recip_rhoConst
275             ELSE             ELSE
276              phiHydC(i,j) = phiHydF(i,j)              phiHydC(i,j) = phiHydF(i,j)
277       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst       &        +dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
# Line 328  C--     Calculate specific volume anomal Line 326  C--     Calculate specific volume anomal
326            DO i=iMin,iMax            DO i=iMin,iMax
327              locAlpha=alphaRho(i,j)+rhoConst              locAlpha=alphaRho(i,j)+rhoConst
328              alphaRho(i,j)=maskC(i,j,k,bi,bj)*              alphaRho(i,j)=maskC(i,j,k,bi,bj)*
329       &              (one/locAlpha - recip_rhoConst)       &              (oneRL/locAlpha - recip_rhoConst)
330            ENDDO            ENDDO
331          ENDDO          ENDDO
332    
333    #ifdef ALLOW_MOM_COMMON
334    C--     Quasi-hydrostatic terms are added as if they modify the specific-volume
335            IF (quasiHydrostatic) THEN
336             CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
337            ENDIF
338    #endif /* ALLOW_MOM_COMMON */
339    
340  C----  Hydrostatic pressure at cell centers  C----  Hydrostatic pressure at cell centers
341    
342         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
# Line 352  C           conserve KE+PE exactly even Line 357  C           conserve KE+PE exactly even
357               phiHydC(i,j) = ddRloc*alphaRho(i,j)               phiHydC(i,j) = ddRloc*alphaRho(i,j)
358  c--to reproduce results of c48d_post: uncomment those 4+1 lines  c--to reproduce results of c48d_post: uncomment those 4+1 lines
359  c            phiHydC(i,j)=phiHydF(i,j)  c            phiHydC(i,j)=phiHydF(i,j)
360  c    &          +(hFacC(i,j,k,bi,bj)-half)*drF(k)*alphaRho(i,j)  c    &          +(hFacC(i,j,k,bi,bj)-halfRL)*drF(k)*alphaRho(i,j)
361  c            phiHydF(i,j)=phiHydF(i,j)  c            phiHydF(i,j)=phiHydF(i,j)
362  c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)  c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
363             ELSE             ELSE
364               phiHydC(i,j) = phiHydF(i,j) + half*drF(k)*alphaRho(i,j)               phiHydC(i,j) = phiHydF(i,j) + halfRL*drF(k)*alphaRho(i,j)
365  c            phiHydF(i,j) = phiHydF(i,j) +      drF(k)*alphaRho(i,j)  c            phiHydF(i,j) = phiHydF(i,j) +        drF(k)*alphaRho(i,j)
366             ENDIF             ENDIF
367  c-- and comment this last one:  c-- and comment this last one:
368               phiHydF(i,j) = phiHydC(i,j) + half*drF(k)*alphaRho(i,j)               phiHydF(i,j) = phiHydC(i,j) + halfRL*drF(k)*alphaRho(i,j)
369  c-----  c-----
370            ENDDO            ENDDO
371           ENDDO           ENDDO
# Line 368  c----- Line 373  c-----
373         ELSE         ELSE
374  C  --  Finite Difference Form, with Part-Cell Bathy  C  --  Finite Difference Form, with Part-Cell Bathy
375    
376           dRlocM=half*drC(k)           dRlocM = halfRL*drC(k)
377           IF (k.EQ.1) dRlocM=rF(k)-rC(k)           IF (k.EQ.1) dRlocM=rF(k)-rC(k)
378           IF (k.EQ.Nr) THEN           IF (k.EQ.Nr) THEN
379             dRlocP=rC(k)-rF(k+1)             dRlocP=rC(k)-rF(k+1)
380           ELSE           ELSE
381             dRlocP=half*drC(k+1)             dRlocP=halfRL*drC(k+1)
382           ENDIF           ENDIF
383           rec_dRm = one/(rF(k)-rC(k))           rec_dRm = oneRL/(rF(k)-rC(k))
384           rec_dRp = one/(rC(k)-rF(k+1))           rec_dRp = oneRL/(rC(k)-rF(k+1))
385    
386           DO j=jMin,jMax           DO j=jMin,jMax
387            DO i=iMin,iMax            DO i=iMin,iMax
# Line 388  C---------- This discretization is the " Line 393  C---------- This discretization is the "
393  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
394               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
395  #endif  #endif
396               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*dRlocM               phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
397       &                      +MIN(zero,ddRloc)*rec_dRp*dRlocP       &                      +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
398       &                     )*alphaRho(i,j)       &                     )*alphaRho(i,j)
399             ELSE             ELSE
400               phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)               phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
# Line 410  C       the specific volume, analogous t Line 415  C       the specific volume, analogous t
415  C--     virtual potential temperature anomaly (including water vapour effect)  C--     virtual potential temperature anomaly (including water vapour effect)
416          DO j=jMin,jMax          DO j=jMin,jMax
417           DO i=iMin,iMax           DO i=iMin,iMax
418            alphaRho(i,j)=maskC(i,j,k,bi,bj)            alphaRho(i,j) = ( tFld(i,j,k,bi,bj)
419       &             *( tFld(i,j,k,bi,bj)*(sFld(i,j,k,bi,bj)*atm_Rq+one)       &                      *( sFld(i,j,k,bi,bj)*atm_Rq + oneRL )
420       &               -tRef(k) )       &                    - tRef(k) )*maskC(i,j,k,bi,bj)
421           ENDDO           ENDDO
422          ENDDO          ENDDO
423    
424    #ifdef ALLOW_MOM_COMMON
425    C--     Quasi-hydrostatic terms are added in as if they modify the Pot.Temp
426            IF (quasiHydrostatic) THEN
427             CALL MOM_QUASIHYDROSTATIC(bi,bj,k,uVel,vVel,alphaRho,myThid)
428            ENDIF
429    #endif /* ALLOW_MOM_COMMON */
430    
431  C---    Integrate d Phi / d pi  C---    Integrate d Phi / d pi
432    
433         IF (integr_GeoPot.EQ.0) THEN         IF (integr_GeoPot.EQ.0) THEN
# Line 431  C--------------------------------------- Line 443  C---------------------------------------
443       &                   -((rC( k )/atm_Po)**atm_kappa) )       &                   -((rC( k )/atm_Po)**atm_kappa) )
444           ELSE           ELSE
445             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
446       &                   -((rC( k )/atm_Po)**atm_kappa) )*half       &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
447           ENDIF           ENDIF
448           IF (k.EQ.Nr) THEN           IF (k.EQ.Nr) THEN
449             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
450       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
451           ELSE           ELSE
452             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
453       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
454           ENDIF           ENDIF
455  C-------- This discretization is the energy conserving form  C-------- This discretization is the energy conserving form
456           DO j=jMin,jMax           DO j=jMin,jMax
# Line 497  C--------- Line 509  C---------
509       &                   -((rC( k )/atm_Po)**atm_kappa) )       &                   -((rC( k )/atm_Po)**atm_kappa) )
510           ELSE           ELSE
511             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)             ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
512       &                   -((rC( k )/atm_Po)**atm_kappa) )*half       &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
513           ENDIF           ENDIF
514           IF (k.EQ.Nr) THEN           IF (k.EQ.Nr) THEN
515             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
516       &                   -((rF(k+1)/atm_Po)**atm_kappa) )       &                   -((rF(k+1)/atm_Po)**atm_kappa) )
517           ELSE           ELSE
518             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)             ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
519       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*half       &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
520           ENDIF           ENDIF
521           rec_dRm = one/(rF(k)-rC(k))           rec_dRm = oneRL/(rF(k)-rC(k))
522           rec_dRp = one/(rC(k)-rF(k+1))           rec_dRp = oneRL/(rC(k)-rF(k+1))
523           DO j=jMin,jMax           DO j=jMin,jMax
524            DO i=iMin,iMax            DO i=iMin,iMax
525             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN             IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
# Line 515  C--------- Line 527  C---------
527  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
528               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)               ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
529  #endif  #endif
530               phiHydC(i,j) =( MAX(zero,ddRloc)*rec_dRm*ddPIm               phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
531       &                      +MIN(zero,ddRloc)*rec_dRp*ddPIp )       &                      +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
532       &                    *alphaRho(i,j)       &                     )*alphaRho(i,j)
533             ELSE             ELSE
534               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)               phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
535             ENDIF             ENDIF
Legend:
Removed from v.1.41  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.22