/[MITgcm]/MITgcm/pkg/opps/opps_calc.F

Diff of /MITgcm/pkg/opps/opps_calc.F

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

-revision 1.10 by mlosch,
Tue Jul 19 13:08:24 2011 UTC
+revision 1.11 by jmc,
Wed Jun 27 22:34:07 2012 UTC
 Line 2 
 C $Header$
  C $Name$
  #include "OPPS_OPTIONS.h"
+ #undef OPPS_ORGCODE
+ C--  File opps_calc.F:
+ C--   Contents
+ C--   o OPPS_CALC
+ C--   o STATE1
+ C--   o NLOPPS
+ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
  CBOP
  C !ROUTINE: OPPS_CALC
-Line 14 
 C !INTERFACE: ==========================
+Line 22 
 C !INTERFACE: ==========================
       I     I, J, bi, bj, myTime, myIter, myThid )
  C !DESCRIPTION: \bv
- C     /=====================================================================\
+ C     *=====================================================================*
  C     | SUBROUTINE OPPS_CALC                                                |
  C     | o Compute all OPPS fields defined in OPPS.h                         |
- C     |=====================================================================|
+ C     *=====================================================================*
  C     | This subroutine is based on the routine 3dconvection.F              |
  C     | by E. Skyllingstad (?)                                              |
  C     | plenty of modifications to make it work:                            |
-Line 53 
 C     |   Td(:) = Pd(:,1), etc.
+Line 61 
 C     |   Td(:) = Pd(:,1), etc.
  C     | o TODO:                                                             |
  C     |   clean up the logic of the vertical loops and get rid off the      |
  C     |   GOTO statements                                                   |
- C     \=====================================================================/
+ C     *=====================================================================*
-       IMPLICIT NONE
- C
- C--------------------------------------------------------------------
  C \ev
+       IMPLICIT NONE
  C !USES: ============================================================
  #include "SIZE.h"
  #include "EEPARAMS.h"
  #include "PARAMS.h"
  #include "OPPS.h"
- #include "FFIELDS.h"
  #include "GRID.h"
-       EXTERNAL DIFFERENT_MULTIPLE
-       LOGICAL  DIFFERENT_MULTIPLE
  C !INPUT PARAMETERS: ===================================================
- c Routine arguments
+ C Routine arguments
- c     bi, bj - array indices on which to apply calculations
+ C     bi, bj :: array indices on which to apply calculations
- c     myTime - Current time in simulation
+ C     myTime :: Current time in simulation
        INTEGER I, J, bi, bj, KMax, nTracer, nTracerInUse
        INTEGER myThid, myIter
-Line 82 
 c     myTime - Current time in simulatio
+Line 84 
 c     myTime - Current time in simulatio
        _RL tracerEnv(Nr,nTracer),wVel(Nr)
  #ifdef ALLOW_OPPS
+ C !FUNCTIONS: ==========================================================
+ c     EXTERNAL DIFFERENT_MULTIPLE
+ c     LOGICAL  DIFFERENT_MULTIPLE
+       _RL STATE1
  C !LOCAL VARIABLES: ====================================================
- c Local constants
+ C Local constants
- C     imin, imax, jmin, jmax  - array computation indices
+ C     msgBuf      :: Informational/error message buffer
- C     msgBuf      - Informational/error message buffer
        CHARACTER*(MAX_LEN_MBUF) msgBuf
        INTEGER K, K2, K2m1, K2p1, ktr
        INTEGER ntime,nn,kmx,ic
-Line 93 
 C     msgBuf      - Informational/error
+Line 99 
 C     msgBuf      - Informational/error
        _RL wsqr,oldflux,newflux,entrainrate
        _RL pmix
-       _RL D1,D2,state1
+       _RL D1, D2
        _RL dz1,dz2
        _RL radius,StartingFlux
        _RL dtts,dt
-Line 101 
 C     Arrays
+Line 107 
 C     Arrays
        _RL Paa(Nr,nTracer)
        _RL wda(Nr), mda(Nr), pda(Nr,nTracer)
  C
- C     Pd, Wd           - tracers, vertical velocity in plume
+ C     Pd, Wd           :: tracers, vertical velocity in plume
- C     Md               - plume mass flux (?)
+ C     Md               :: plume mass flux (?)
- C     Ad               - fractional area covered by plume
+ C     Ad               :: fractional area covered by plume
- C     Dd               - density in plume
+ C     Dd               :: density in plume
- C     De               - density of environment
+ C     De               :: density of environment
- C     PlumeEntrainment -
+ C     PlumeEntrainment ::
        _RL Ad(Nr),Wd(Nr),Dd(Nr),Md(Nr)
        _RL De(Nr)
        _RL PlumeEntrainment(Nr)
        _RL Pd(Nr,nTracer)
  CEOP
  C--   Check to see if should convect now
  C      IF ( DIFFERENT_MULTIPLE(cAdjFreq,myTime,deltaTClock) ) THEN
        IF ( .true. ) THEN
-Line 121 
 C     local initialization
+Line 126 
 C     local initialization
  C     Copy some arrays
        dtts = dTtracerLev(1)
- C
  C     start k-loop
- C
        DO k=1,KMax-1
- c
- c initialize the plume T,S,density, and w velocity
+ C initialize the plume T,S,density, and w velocity
- c
         DO ktr=1,nTracerInUse
          Pd(k,ktr) = tracerEnv(k,ktr)
         ENDDO
-        Dd(k)=state1(Pd(k,2),Pd(k,1),i,j,k,bi,bj,myThid)
+        Dd(k)=STATE1(Pd(k,2),Pd(k,1),i,j,k,bi,bj,myThid)
         De(k)=Dd(k)
  CML       print *, 'ml-opps:', i,j,k,tracerEnv(k,2),tracerEnv(k,1),
  CML     &      Dd(k),Pd(k,1),Pd(k,2)
-Line 142 
 CML(
+Line 146 
 CML(
  CML    avoid division by zero
  CML       IF (Wd(K) .EQ. 0.D0) Wd(K) = 2.23e-16
  CML)
- c
- c guess at initial top grid cell vertical velocity
+ C guess at initial top grid cell vertical velocity
- c
  CML          Wd(k) = 0.03
- c
- c these estimates of initial plume velocity based on plume size and
+ C these estimates of initial plume velocity based on plume size and
- c top grid cell water mass
+ C top grid cell water mass
- c
  c          Wd(k) = 0.5*drF(k)/(dtts*FRACTIONAL_AREA)
  c          Wd(k) = 0.5*drF(k)/dtts
- c
         wsqr=Wd(k)*Wd(k)
         PlumeEntrainment(k) = 0.0
- c
- c
- c
  #ifdef ALLOW_OPPS_DEBUG
         IF ( OPPSdebugLevel.GE.debLevB ) THEN
          WRITE(msgBuf,'(A,I3)')
-Line 172 
 c
+Line 174 
 c
         dz2=DrF(k)
         DO k2=k,KMax-1
-         D1=state1( Pd(k2,2), Pd(k2,1),i,j,k2+1,bi,bj,myThid)
+         D1=STATE1( Pd(k2,2), Pd(k2,1),i,j,k2+1,bi,bj,myThid)
-         D2=state1( tracerEnv(k2+1,2), tracerEnv(k2+1,1),
+         D2=STATE1( tracerEnv(k2+1,2), tracerEnv(k2+1,1),
       &                                i,j,k2+1,bi,bj,myThid)
          De(k2+1)=D2
- c
- c To start downward, parcel has to initially be heavier than environment
+ C To start downward, parcel has to initially be heavier than environment
- c but after it has started moving, we continue plume until plume tke or
+ C but after it has started moving, we continue plume until plume tke or
- c flux goes negative
+ C flux goes negative
- c
  CML     &       _hFacC(i,j,k-1,bi,bj)
  CML     &       *_hFacC(i,j,k,bi,bj) .GT. 0.
  CML     &  .AND.
          IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN
           dz1=dz2
           dz2=DrF(k2+1)
- c
  C     find mass flux according to eq.(3) from paper by vertical integration
- c
           newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*
       &        .5*(dz1+dz2)
  CML         print *, 'ml-opps:', i,j,k,oldflux,newflux,e2,radius,
  CML     &        Wd(k2),Dd(k2),Pd(k2,1),Pd(k2,2),dz1,dz2
- c
           PlumeEntrainment(k2+1) = newflux/StartingFlux
- c
           IF(newflux.LE.0.0) then
  #ifdef ALLOW_OPPS_DEBUG
            IF ( OPPSdebugLevel.GE.debLevA ) THEN
-Line 210 
 c
+Line 212 
 c
            if(maxdepth.eq.k) goto 1000
            goto 1
           endif
- c
- c entrainment rate is basically a scaled mass flux dM/M
+ C entrainment rate is basically a scaled mass flux dM/M
- c
           entrainrate = (newflux - oldflux)/newflux
           oldflux = newflux
- c
- c
+ C mix var(s) are the average environmental values over the two grid levels
- c mix var(s) are the average environmental values over the two grid levels
- c
           DO ktr=1,nTracerInUse
            pmix=(dz1*tracerEnv(k2,ktr)+dz2*tracerEnv(k2+1,ktr))
       &         /(dz1+dz2)
            Pd(k2+1,ktr)=Pd(k2,ktr)
       &         - entrainrate*(pmix - Pd(k2,ktr))
           ENDDO
- c
- c compute the density at this level for the buoyancy term in the
+ C compute the density at this level for the buoyancy term in the
- c vertical k.e. equation
+ C vertical k.e. equation
- c
-          Dd(k2+1)=state1(Pd(k2+1,2),Pd(k2+1,1),i,j,k2+1,bi,bj,myThid)
+          Dd(k2+1)=STATE1(Pd(k2+1,2),Pd(k2+1,1),i,j,k2+1,bi,bj,myThid)
- c
- c next, solve for the vertical velocity k.e. using combined eq. (4)
+ C next, solve for the vertical velocity k.e. using combined eq. (4)
- c and eq (5) from the paper
+ C and eq (5) from the paper
- c
  #ifdef ALLOW_OPPS_DEBUG
           IF ( OPPSdebugLevel.GE.debLevA ) THEN
            WRITE(msgBuf,'(A,3E12.4,I3)')
-Line 246 
 CML   insert Eq. (4) into Eq. (5) to get
+Line 247 
 CML   insert Eq. (4) into Eq. (5) to get
           wsqr = wsqr - wsqr*abs(entrainrate)+ gravity*
       &        (dz1*(Dd(k2)-De(k2))/De(k2)
       &        +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1))
- c
- c if negative k.e. then plume has reached max depth, get out of loop
+ C if negative k.e. then plume has reached max depth, get out of loop
- c
           IF(wsqr.LE.0.0)then
            maxdepth = k2
  #ifdef ALLOW_OPPS_DEBUG
-Line 273 
 c
+Line 274 
 c
            goto 1
           endif
           Wd(k2+1)=sqrt(wsqr)
- C
  C     compute a new radius based on the new mass flux at this grid level
  C     from Eq. (4)
- C
           radius=sqrt(newflux/(Wd(k2)*Dd(k2)))
          ELSE
           maxdepth=k2
-Line 284 
 C
+Line 285 
 C
           GOTO 1
          ENDIF
         ENDDO
- c
- c plume has reached the bottom
+ C plume has reached the bottom
- c
         MaxDepth=kMax
- c
     CONTINUE
- c
         Ad(k)=FRACTIONAL_AREA
         IC=0
- c
- c start iteration on fractional area, not used in OGCM implementation
+ C start iteration on fractional area, not used in OGCM implementation
- c
- c
         DO IC=1,Max_ABE_Iterations
- c
- c
+ C next compute the mass flux beteen each grid box using the entrainment
- c next compute the mass flux beteen each grid box using the entrainment
- c
          Md(k)=Wd(k)*Ad(k)
- c
          DO k2=k+1,maxDepth
           Md(k2)=Md(k)*PlumeEntrainment(k2)
  #ifdef ALLOW_OPPS_DEBUG
-Line 315 
 c
+Line 314 
 c
           ENDIF
  #endif /* ALLOW_OPPS_DEBUG */
          ENDDO
- c
- c Now move on to calculate new temperature using flux from
+ C Now move on to calculate new temperature using flux from
- c Td, Sd, Wd, ta, sa, and we. Values for these variables are at
+ C Td, Sd, Wd, ta, sa, and we. Values for these variables are at
- c center of grid cell, use weighted average to get boundary values
+ C center of grid cell, use weighted average to get boundary values
- c
- c use a timestep limited by the GCM model timestep and the maximum plume
+ C use a timestep limited by the GCM model timestep and the maximum plume
- c velocity (CFL criteria)
+ C velocity (CFL criteria)
- c
- c
+ C calculate the weighted wd, td, and sd
- c calculate the weighted wd, td, and sd
- c
          dt = dtts
          do k2=k,maxDepth-1
           IF ( Wd(K2) .NE. 0. _d 0 ) dt = min(dt,drF(k2)/Wd(k2))
- c
- c time integration will be integer number of steps to get one
+ C time integration will be integer number of steps to get one
- c gcm time step
+ C gcm time step
- c
           ntime = nint(0.5*int(dtts/dt))
           if(ntime.eq.0) then
            ntime = 1
           endif
- c
- c make sure area weighted vertical velocities match; in other words
+ C make sure area weighted vertical velocities match; in other words
- c make sure mass in equals mass out at the intersection of each grid
+ C make sure mass in equals mass out at the intersection of each grid
- c cell. Eq. (20)
+ C cell. Eq. (20)
- c
           mda(k2) = (md(k2)*drF(k2)+md(k2+1)*drF(k2+1))/
       *        (drF(k2)+drF(k2+1))
- c
           wda(k2) = (wd(k2)*drF(k2)+wd(k2+1)*drF(k2+1))/
       *        (drF(k2)+drF(k2+1))
- c
           DO ktr = 1, nTracerInUse
            Pda(k2,ktr) = Pd(k2,ktr)
            Paa(k2,ktr) = tracerEnv(k2+1,ktr)
           ENDDO
- c
          enddo
          dt = min(dt,dtts)
  #ifdef ALLOW_OPPS_DEBUG
-Line 366 
 c
+Line 363 
 c
          DO ktr=1,nTracerInUse
           Pda(maxdepth,ktr) = Pd(maxdepth,ktr)
          ENDDO
- C
          kmx = maxdepth-1
          do nn=1,ntime
- C
  C     top point
- C
           DO ktr = 1,nTracerInUse
            tracerEnv(k,ktr) =  tracerEnv(k,ktr)-
       &        (mda(k)*(Pda(k,ktr)-Paa(k,ktr)))*dt*recip_drF(k)
           ENDDO
- c
- c now do inner points if there are any
+ C now do inner points if there are any
- c
  CML         if(Maxdepth-k.gt.1) then
  CML    This if statement is superfluous
  CML         IF ( k .LT. Maxdepth-1 ) THEN
-Line 387 
 CML         mda(maxDepth) = 0.
+Line 384 
 CML         mda(maxDepth) = 0.
           DO k2=k+1,kmx
            k2m1 = max(k,k2-1)
            k2p1 = max(k2+1,maxDepth)
- c
             DO ktr = 1,nTracerInUse
              tracerEnv(k2,ktr) = tracerEnv(k2,ktr) +
       &           (mda(k2m1)*(Pda(k2m1,ktr)-Paa(k2m1,ktr))
-Line 397 
 c
+Line 394 
 c
            ENDDO
  CML    This if statement is superfluous
  CML         ENDIF
- C
  C     bottom point
- C
           DO ktr=1,nTracerInUse
            tracerEnv(kmx+1,ktr) =  tracerEnv(kmx+1,ktr)+
       &        mda(kmx)*(Pda(kmx,ktr)-Paa(kmx,ktr))*dt*recip_drF(kmx+1)
           ENDDO
- c
- c     set the environmental temp and salinity to equal new fields
+ C     set the environmental temp and salinity to equal new fields
- c
           DO ktr=1,nTracerInUse
            DO k2=1,kmx
             paa(k2,ktr) = tracerEnv(k2+1,ktr)
            ENDDO
           ENDDO
- c
- c end loop on number of time integration steps
+ C end loop on number of time integration steps
- c
          enddo
         ENDDO
   continue
- C
  C     count convection event in this grid cell
- C
         OPPSconvectCount(I,J,K,bi,bj) =
       &      OPPSconvectCount(I,J,K,bi,bj) + 1. _d 0
- C
  C     jump here if k = maxdepth or if level not unstable, go to next
  C     profile point
- C
  continue
- c
- C
  C     end  of k-loop
- C
        ENDDO
  C--   End IF (DIFFERENT_MULTIPLE)
-Line 439 
 C--   End IF (DIFFERENT_MULTIPLE)
+Line 435 
 C--   End IF (DIFFERENT_MULTIPLE)
        RETURN
        END
-       _RL FUNCTION STATE1(sLoc,tLoc,I,J,KREF,bi,bj,mythid)
+ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
+       _RL FUNCTION STATE1( sLoc, tLoc, i, j, kRef, bi, bj, myThid )
  C     !DESCRIPTION: \bv
  C     *===============================================================*
  C     | o SUBROUTINE STATE1
 Line 456 
 C     == Global variables ==
  #include "SIZE.h"
  #include "EEPARAMS.h"
  #include "PARAMS.h"
- #include "EOS.h"
  #include "GRID.h"
  #include "DYNVARS.h"
  C     !INPUT/OUTPUT PARAMETERS:
  C     == Routine arguments ==
-       INTEGER I,J,kRef,bi,bj,myThid
+       INTEGER i, j, kRef, bi, bj, myThid
-       _RL tLoc,sLoc
+       _RL tLoc, sLoc
  C     !LOCAL VARIABLES:
  C     == Local variables ==
-       _RL rhoLoc, dRho
+       _RL rhoLoc
        _RL pLoc
-       _RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1
-       _RL ct, sa, sqrtsa, p
-       _RL rfresh, rsalt, rhoP0
-       _RL bMfresh, bMsalt, bMpres, BulkMod
-       _RL rhoNum, rhoDen, den, epsln
-       PARAMETER ( epsln = 0.D0 )
-       character*(max_len_mbuf) msgbuf
  CMLC     estimate pressure from depth at cell centers
  CML      mtoSI = gravity*rhoConst
  CML      pLoc = ABS(rC(kRef))*mtoSI
-       IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
+       IF ( usingZCoords ) THEN
  C     in Z coordinates the pressure is rho0 * (hydrostatic) Potential
         IF ( useDynP_inEos_Zc ) THEN
  C----------
-Line 496 
 C----------
+Line 487 
 C----------
         ELSE
          pLoc = -rhoConst*rC(kRef)*gravity*maskC(i,j,kRef,bi,bj)
         ENDIF
-       ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
+       ELSEIF ( usingPCoords ) THEN
  C     in P coordinates the pressure is just the coordinate of
  C     the tracer point
         pLoc = rC(kRef)* maskC(i,j,kRef,bi,bj)
        ENDIF
-       rhoLoc  = 0. _d 0
+       CALL FIND_RHO_SCALAR( tLoc, sLoc, pLoc, rhoLoc, myThid )
-       rhoP0   = 0. _d 0
+       STATE1 = rhoLoc
-       bulkMod = 0. _d 0
-       rfresh  = 0. _d 0
-       rsalt   = 0. _d 0
-       bMfresh = 0. _d 0
-       bMsalt  = 0. _d 0
-       bMpres  = 0. _d 0
-       rhoNum  = 0. _d 0
-       rhoDen  = 0. _d 0
-       den     = 0. _d 0
-       t1 = tLoc
-       t2 = t1*t1
-       t3 = t2*t1
-       t4 = t3*t1
-       s1  = sLoc
-       IF ( equationOfState .EQ. 'LINEAR' ) THEN
-        dRho = rhoNil-rhoConst
-        rhoLoc=rhoNil* (
-      &        sBeta *(sLoc-sRef(kRef))
-      &      - tAlpha*(tLoc-tRef(KREF)) ) + dRho
-       ELSEIF (equationOfState.EQ.'POLY3') THEN
- C     this is not correct, there is a field eosSig0 which should be use here
- C     but I DO not intent to include the reference level in this routine
-          WRITE(*,'(a)')
-      &        ' FIND_RHO_SCALAR: for POLY3, the density is not'
-          WRITE(*,'(a)')
-      &         '                 computed correctly in this routine'
-          rhoLoc = 0. _d 0
-       ELSEIF ( equationOfState(1:5).EQ.'JMD95'
-      &      .OR. equationOfState.EQ.'UNESCO' ) THEN
- C     nonlinear equation of state in pressure coordinates
-          s3o2 = s1*SQRT(s1)
-          p1 = pLoc*SItoBar
-          p2 = p1*p1
- C     density of freshwater at the surface
-          rfresh =
-      &          eosJMDCFw(1)
-      &        + eosJMDCFw(2)*t1
-      &        + eosJMDCFw(3)*t2
-      &        + eosJMDCFw(4)*t3
-      &        + eosJMDCFw(5)*t4
-      &        + eosJMDCFw(6)*t4*t1
- C     density of sea water at the surface
-          rsalt =
-      &        s1*(
-      &             eosJMDCSw(1)
-      &           + eosJMDCSw(2)*t1
-      &           + eosJMDCSw(3)*t2
-      &           + eosJMDCSw(4)*t3
-      &           + eosJMDCSw(5)*t4
-      &           )
-      &        + s3o2*(
-      &             eosJMDCSw(6)
-      &           + eosJMDCSw(7)*t1
-      &           + eosJMDCSw(8)*t2
-      &           )
-      &           + eosJMDCSw(9)*s1*s1
-          rhoP0 = rfresh + rsalt
- C     secant bulk modulus of fresh water at the surface
-          bMfresh =
-      &             eosJMDCKFw(1)
-      &           + eosJMDCKFw(2)*t1
-      &           + eosJMDCKFw(3)*t2
-      &           + eosJMDCKFw(4)*t3
-      &           + eosJMDCKFw(5)*t4
- C     secant bulk modulus of sea water at the surface
-          bMsalt =
-      &        s1*( eosJMDCKSw(1)
-      &           + eosJMDCKSw(2)*t1
-      &           + eosJMDCKSw(3)*t2
-      &           + eosJMDCKSw(4)*t3
-      &           )
-      &    + s3o2*( eosJMDCKSw(5)
-      &           + eosJMDCKSw(6)*t1
-      &           + eosJMDCKSw(7)*t2
-      &           )
- C     secant bulk modulus of sea water at pressure p
-          bMpres =
-      &        p1*( eosJMDCKP(1)
-      &           + eosJMDCKP(2)*t1
-      &           + eosJMDCKP(3)*t2
-      &           + eosJMDCKP(4)*t3
-      &           )
-      &   + p1*s1*( eosJMDCKP(5)
-      &           + eosJMDCKP(6)*t1
-      &           + eosJMDCKP(7)*t2
-      &           )
-      &      + p1*s3o2*eosJMDCKP(8)
-      &      + p2*( eosJMDCKP(9)
-      &           + eosJMDCKP(10)*t1
-      &           + eosJMDCKP(11)*t2
-      &           )
-      &    + p2*s1*( eosJMDCKP(12)
-      &           + eosJMDCKP(13)*t1
-      &           + eosJMDCKP(14)*t2
-      &           )
-          bulkMod = bMfresh + bMsalt + bMpres
- C     density of sea water at pressure p
-          rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst
-       ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
-          sp5 = SQRT(s1)
-          p1   = pLoc*SItodBar
-          p1t1 = p1*t1
-          rhoNum = eosMDJWFnum(0)
-      &        + t1*(eosMDJWFnum(1)
-      &        +     t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
-      &        + s1*(eosMDJWFnum(4)
-      &        +     eosMDJWFnum(5)*t1  + eosMDJWFnum(6)*s1)
-      &        + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
-      &        +     eosMDJWFnum(9)*s1
-      &        +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
-          den = eosMDJWFden(0)
-      &        + t1*(eosMDJWFden(1)
-      &        +     t1*(eosMDJWFden(2)
-      &        +         t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
-      &        + s1*(eosMDJWFden(5)
-      &        +     t1*(eosMDJWFden(6)
-      &        +         eosMDJWFden(7)*t2)
-      &        +     sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
-      &        + p1*(eosMDJWFden(10)
-      &        +     p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
-          rhoDen = 1.0/(epsln+den)
-          rhoLoc = rhoNum*rhoDen - rhoConst
-       ELSEIF( equationOfState .EQ. 'TEOS10' ) THEN
-        ct      = tLoc
-        sa      = sLoc
-        IF ( sa .GT. 0. _d 0 ) THEN
-         sqrtsa = SQRT(sa)
-        ELSE
-         sa     = 0. _d 0
-         sqrtsa = 0. _d 0
-        ENDIF
-        p       = pLoc*SItodBar
-        rhoNum = teos(01)
-      &   + ct*(teos(02) + ct*(teos(03) + teos(04)*ct))
-      &   + sa*(teos(05) + ct*(teos(06) + teos(07)*ct)
-      &   + sqrtsa*(teos(08) + ct*(teos(09)
-      &            + ct*(teos(10) + teos(11)*ct))))
-      &   + p*(teos(12) + ct*(teos(13) + teos(14)*ct)
-      &                      + sa*(teos(15) + teos(16)*ct)
-      &   + p*(teos(17) + ct*(teos(18) + teos(19)*ct) + teos(20)*sa))
-        den = teos(21)
-      &   + ct*(teos(22) + ct*(teos(23) + ct*(teos(24) + teos(25)*ct)))
-      &   + sa*(teos(26) + ct*(teos(27) + ct*(teos(28)
-      &   + ct*(teos(29) + teos(30)*ct)))
-      &   + teos(36)*sa
-      %   + sqrtsa*(teos(31) + ct*(teos(32) + ct*(teos(33)
-      &            + ct*(teos(34) + teos(35)*ct)))))
-      %   + p*(teos(37) + ct*(teos(38) + ct*(teos(39) + teos(40)*ct))
-      %   + sa*(teos(41) + teos(42)*ct)
-      %   + p*(teos(43) + ct*(teos(44) + teos(45)*ct + teos(46)*sa)
-      %   + p*(teos(47) + teos(48)*ct)))
-        rhoDen = 1.0/(epsln+den)
-        rhoLoc = rhoNum*rhoDen
-       ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
- C
-       ELSE
-        WRITE(msgbuf,'(3A)')
-      &        ' STATE1 : equationOfState = "',
-      &        equationOfState,'"'
-        CALL PRINT_ERROR( msgbuf, mythid )
-        STOP 'ABNORMAL END: S/R STATE1 in OPPS_CALC'
-       ENDIF
-       state1 = rhoLoc + rhoConst
  #endif /* ALLOW_OPPS */
        RETURN
        END
- #undef OPPS_ORGCODE
  #ifdef OPPS_ORGCODE
- c Listed below is the subroutine for use in parallel 3-d circulation code.
+ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
- c It has been used in the parallel semtner-chervin code and is now being used
- c In the POP code.  The subroutine is called nlopps (long story to explain why).
+ C Listed below is the subroutine for use in parallel 3-d circulation code.
+ C It has been used in the parallel semtner-chervin code and is now being used
- c I've attached the version of lopps that we've been using in the simulations.
+ C In the POP code.  The subroutine is called nlopps (long story to explain why).
- c There is one common block that is different from the standard model commons
- c (countc) and it is not needed if the array convadj is not used.  The routine
+ C I've attached the version of lopps that we've been using in the simulations.
- c does need "kmp" which is why the boundc common is included. For massively
+ C There is one common block that is different from the standard model commons
- c parallel codes (like POP) we think this will work well when converted from a
+ C (countc) and it is not needed if the array convadj is not used.  The routine
- c "slab" (i=is,ie) to a column, which just means removing the "do i=is,ie" loop.  c There are differences between this
+ C does need "kmp" which is why the boundc common is included. For massively
- c code and the 1-d code and the earlier scheme implemented in 3-d models. These c differences are described below.
+ C parallel codes (like POP) we think this will work well when converted from a
+ C "slab" (i=is,ie) to a column, which just means removing the "do i=is,ie" loop.  c There are differences between this
+ C code and the 1-d code and the earlier scheme implemented in 3-d models. These c differences are described below.
        subroutine nlopps(j,is,ie,ta,sa,gcmdz)
- c
        parameter (imt = 361 , jmt = 301 , km = 30 )
- c
- c     Nlopps:   E. Skyllingstad and T. Paluszkiewicz
- c
- c     Version: December 11, 1996
- c
- c     Nlopps:  This version of lopps is significantly different from
- c     the original code developed by R. Romea and T. Paluskiewicz.  The
- c     code uses a flux constraint to control the change in T and S at
- c     each grid level.  First, a plume profile of T,S, and W are
- c     determined using the standard plume model, but with a detraining
- c     mass instead of entraining.  Thus, the T and S plume
- c     characteristics still change, but the plume contracts in size
- c     rather than expanding ala classical entraining plumes.  This
- c     is heuristically more in line with large eddy simulation results.
- c     At each grid level, the convergence of plume velocity determines
- c     the flux of T and S, which is conserved by using an upstream
- c     advection.  The vertical velocity is balanced so that the area
- c     weighted upward velocity equals the area weighted downdraft
- c     velocity, ensuring mass conservation. The present implementation
- c     adjusts the plume for a time period equal to the time for 1/2 of
- c     the mass of the fastest moving level to move downward.  As a
- c     consequence, the model does not completely adjust the profile at
- c     each model time step, but provides a smooth adjustment over time.
- c
- c
+ C     Nlopps:   E. Skyllingstad and T. Paluszkiewicz
- c
+ C     Version: December 11, 1996
+ C     Nlopps:  This version of lopps is significantly different from
+ C     the original code developed by R. Romea and T. Paluskiewicz.  The
+ C     code uses a flux constraint to control the change in T and S at
+ C     each grid level.  First, a plume profile of T,S, and W are
+ C     determined using the standard plume model, but with a detraining
+ C     mass instead of entraining.  Thus, the T and S plume
+ C     characteristics still change, but the plume contracts in size
+ C     rather than expanding ala classical entraining plumes.  This
+ C     is heuristically more in line with large eddy simulation results.
+ C     At each grid level, the convergence of plume velocity determines
+ C     the flux of T and S, which is conserved by using an upstream
+ C     advection.  The vertical velocity is balanced so that the area
+ C     weighted upward velocity equals the area weighted downdraft
+ C     velocity, ensuring mass conservation. The present implementation
+ C     adjusts the plume for a time period equal to the time for 1/2 of
+ C     the mass of the fastest moving level to move downward.  As a
+ C     consequence, the model does not completely adjust the profile at
+ C     each model time step, but provides a smooth adjustment over time.
  c      include "params.h"
  c      include "plume_fast_inc.h"
-Line 759 
 c #include "loppsd.h"
+Line 551 
 c #include "loppsd.h"
        real ta(imt,km),sa(imt,km),gcmdz(km),dz(km)
        real pdensity,wsqr,oldflux,newflux,entrainrate,adtemp
        REAL Del,D,dza1,dza2,kd,kd1,Smix,Thmix,PlumeS,PlumeT,PlumeD
- c
- c
        INTEGER i,j,k
- clfh
+ Clfh
        integer is,ie,k2
- clfh
+ Clfh
        REAL D1,D2,state1,Density
        REAL dz1,dz2
        REAL radius,StartingFlux
-Line 773 
 clfh
+Line 563 
 clfh
        real wda(km),tda(km),sda(km),mda(km)
        real dtts,dt,sumo,sumn
        integer ntime,nn,kmx,ic
- c
- c
        LOGICAL debug,done
        INTEGER MAX_ABE_ITERATIONS
        PARAMETER(MAX_ABE_ITERATIONS=1)
-Line 798 
 c
+Line 587 
 c
        REAL PlumeEntrainment(km)
        REAL GridThickness(km)
- c
+ C input kmp through a common block
- c input kmp through a common block
- c
        common / boundc / wsx(imt,jmt),wsy(imt,jmt),hfs(imt,jmt),
                  ple(imt,jmt),kmp(imt,jmt),kmq(imt,jmt)
  cwmseas
       &                 ,wsx1(imt,jmt),wsy1(imt,jmt)
                 ,wsx2(imt,jmt),wsy2(imt,jmt)
- c
- c input the variables through a common
+ C input the variables through a common
- c
        logical problem
        common /countc/ convadj(imt,jmt,km),ics,depth(km),problem
- c-----may want to setup an option to get this only on first call
+ C-----may want to setup an option to get this only on first call
- c     otherwise it is repetive
+ C     otherwise it is repetive
- c     griddz is initialize by call to setupgrid
+ C     griddz is initialize by call to setupgrid
- c
- c
          dtts = 2400
- c
          do k=1,km
            dz(k) = 0.01*gcmdz(k)
          enddo
- c
          do k=1,km
             GridThickness(k) = dz(k)
          enddo
- c
- c modified to loop over slab
+ C modified to loop over slab
- c
        DO i=is,ie
          numgridpoints=kmp(i,j)
- c
- c  go to next column if only 1 grid point or on land
+ C  go to next column if only 1 grid point or on land
- c
          if(numgridpoints.le.1) goto 1100
- c
- c loop over depth
+ C loop over depth
- c
- c
        debug = .false.
- c
- c first save copy of initial profile
+ C first save copy of initial profile
- c
        DO k=1,NumGridPoints
           stemp(k)=sa(i,k)
           ttemp(k)=ta(i,k)
- c
- c do a check of t and s range, if out of bounds set flag
+ C do a check of t and s range, if out of bounds set flag
- c
           if(problem) then
              write(0,*)"Code in trouble before this nlopps call"
              return
           endif
- c
           if(sa(i,k).gt.40..or.ta(i,k).lt.-4.0) then
              problem = .true.
              write(0,*)"t out of range at j ",j
-Line 871 
 c
+Line 656 
 c
        endif
        DO k=1,NumGridPoints-1
- c
- c initialize the plume T,S,density, and w velocity
+ C initialize the plume T,S,density, and w velocity
- c
            Sd(k)=stemp(k)
            Td(k)=ttemp(k)
            Dd(k)=state1(stemp(k),ttemp(k),k)
            De(k)=Dd(k)
  c          Wd(k)=VERTICAL_VELOCITY
- c
- c guess at initial top grid cell vertical velocity
+ C guess at initial top grid cell vertical velocity
- c
            Wd(k) = 0.03
- c
- c these estimates of initial plume velocity based on plume size and
+ C these estimates of initial plume velocity based on plume size and
- c top grid cell water mass
+ C top grid cell water mass
- c
  c          Wd(k) = 0.5*dz(k)/(dtts*FRACTIONAL_AREA)
  c          Wd(k) = 0.5*dz(k)/dtts
- c
            wsqr=Wd(k)*Wd(k)
            PlumeEntrainment(k) = 0.0
- c
- c
- c
            if(debug) write(0,*) 'Doing old lowerparcel'
            radius=PlumeRadius
            StartingFlux=radius*radius*Wd(k)*Dd(k)
-Line 905 
 c
+Line 688 
 c
              D1=state1(Sd(k2),Td(k2),k2+1)
              D2=state1(stemp(k2+1),ttemp(k2+1),k2+1)
              De(k2+1)=D2
- c
- c To start downward, parcel has to initially be heavier than environment
+ C To start downward, parcel has to initially be heavier than environment
- c but after it has started moving, we continue plume until plume tke or
+ C but after it has started moving, we continue plume until plume tke or
- c flux goes negative
+ C flux goes negative
- c
              IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN
                   dz1=dz2
                   dz2=GridThickness(k2+1)
- c
- c define mass flux according to eq. 4 from paper
+ C define mass flux according to eq. 4 from paper
- c
                   newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*0.50*
       .              (dz1+dz2)
- c
                   PlumeEntrainment(k2+1) = newflux/StartingFlux
- c
                   IF(newflux.LT.0.0) then
                       if(debug) then
                        write(0,*)"Plume entrained to zero at ",k2
-Line 929 
 c
+Line 712 
 c
                       if(maxdepth.eq.k) goto 1000
                       goto 1
                   endif
- c
- c entrainment rate is basically a scaled mass flux dM/M
+ C entrainment rate is basically a scaled mass flux dM/M
- c
                   entrainrate = (newflux - oldflux)/newflux
                   oldflux = newflux
- c
- c
+ C mix var(s) are the average environmental values over the two grid levels
- c mix var(s) are the average environmental values over the two grid levels
- c
                   smix=(dz1*stemp(k2)+dz2*stemp(k2+1))/(dz1+dz2)
                   thmix=(dz1*ttemp(k2)+dz2*ttemp(k2+1))/(dz1+dz2)
- c
- c first compute the new salinity and temperature for this level
+ C first compute the new salinity and temperature for this level
- c using equations 3.6 and 3.7 from the paper
+ C using equations 3.6 and 3.7 from the paper
- c
- c
- c
                    sd(k2+1)=sd(k2) - entrainrate*(smix - sd(k2))
                    td(k2+1)=td(k2) - entrainrate*(thmix - td(k2))
- c
- c
+ C compute the density at this level for the buoyancy term in the
- c compute the density at this level for the buoyancy term in the
+ C vertical k.e. equation
- c vertical k.e. equation
- c
                   Dd(k2+1)=state1(Sd(k2+1),Td(k2+1),k2+1)
- c
- c next, solve for the vertical velocity k.e. using combined eq. 4
+ C next, solve for the vertical velocity k.e. using combined eq. 4
- c and eq 5 from the paper
+ C and eq 5 from the paper
- c
                   if(debug) then
                    write(0,*)"Dd,De,entr,k ",Dd(k2),De(k2),entrainrate,k2
                   endif
- c
                   wsqr = wsqr - wsqr*abs(entrainrate)+ 9.81*
       .             (dz1*(Dd(k2)-De(k2))/De(k2)
       .             +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1))
- c
- c if negative k.e. then plume has reached max depth, get out of loop
+ C if negative k.e. then plume has reached max depth, get out of loop
- c
                   IF(wsqr.LT.0.0)then
                       maxdepth = k2
                       if(debug) then
-Line 977 
 c
+Line 756 
 c
                       goto 1
                   endif
                   Wd(k2+1)=sqrt(wsqr)
- c
- c compute a new radius based on the new mass flux at this grid level
+ C compute a new radius based on the new mass flux at this grid level
- c
                   radius=sqrt(newflux/(Wd(k2)*Dd(k2)))
                ELSE
                   maxdepth=k2
-Line 987 
 c
+Line 766 
 c
                   GOTO 1
                ENDIF
            ENDDO
- c
- c plume has reached the bottom
+ C plume has reached the bottom
- c
            MaxDepth=NumGridPoints
- c
         continue
- c
            Ad(k)=FRACTIONAL_AREA
            IC=0
- c
- c start iteration on fractional area, not used in OGCM implementation
+ C start iteration on fractional area, not used in OGCM implementation
- c
- c
            DO IC=1,Max_ABE_Iterations
- c
- c
+ C next compute the mass flux beteen each grid box using the entrainment
- c next compute the mass flux beteen each grid box using the entrainment
- c
          continue
               Md(k)=Wd(k)*Ad(k)
- c
               DO k2=k+1,MaxDepth
                 Md(k2)=Md(k)*PlumeEntrainment(k2)
                 if(debug) then
                   write(0,*)"Md, Wd, and  k are ",Md(k2),Wd(k2),k2
                 endif
               ENDDO
- c
- c Now move on to calculate new temperature using flux from
+ C Now move on to calculate new temperature using flux from
- c Td, Sd, Wd, ta, sa, and we. Values for these variables are at
+ C Td, Sd, Wd, ta, sa, and we. Values for these variables are at
- c center of grid cell, use weighted average to get boundary values
+ C center of grid cell, use weighted average to get boundary values
- c
- c use a timestep limited by the GCM model timestep and the maximum plume
+ C use a timestep limited by the GCM model timestep and the maximum plume
- c velocity (CFL criteria)
+ C velocity (CFL criteria)
- c
- c
+ C calculate the weighted wd, td, and sd
- c calculate the weighted wd, td, and sd
- c
               dt = dtts
               do k2=k,maxdepth-1
                  dt = min(dt,dz(k2)/wd(k2))
- c
- c time integration will be integer number of steps to get one
+ C time integration will be integer number of steps to get one
- c gcm time step
+ C gcm time step
- c
                  ntime = nint(0.5*int(dtts/dt))
                  if(ntime.eq.0) then
                     ntime = 1
                  endif
- c
- c make sure area weighted vertical velocities match; in other words
+ C make sure area weighted vertical velocities match; in other words
- c make sure mass in equals mass out at the intersection of each grid
+ C make sure mass in equals mass out at the intersection of each grid
- c cell.
+ C cell.
- c
                  mda(k2) = (md(k2)*dz(k2)+md(k2+1)*dz(k2+1))/
       *                    (dz(k2)+dz(k2+1))
- c
                  wda(k2) = (wd(k2)*dz(k2)+wd(k2+1)*dz(k2+1))/
       *                    (dz(k2)+dz(k2+1))
- c
                  tda(k2) = td(k2)
                  sda(k2) = sd(k2)
- c
                  taa(k2) = ttemp(k2+1)
                  saa(k2) = stemp(k2+1)
- c
               enddo
               dt = min(dt,dtts)
               if(debug) then
-Line 1060 
 c
+Line 836 
 c
               endif
               tda(maxdepth) = td(maxdepth)
               sda(maxdepth) = sd(maxdepth)
- c
- c do top and bottom points first
+ C do top and bottom points first
- c
               kmx = maxdepth-1
               do nn=1,ntime
                 ttemp(k) =  ttemp(k)-
       *                  (mda(k)*(tda(k)-taa(k)))*dt*recip_drF(k)
- c
                 stemp(k) =  stemp(k)-
       *                  (mda(k)*(sda(k)-saa(k)))*dt*recip_drF(k)
- c
- c
+ C now do inner points if there are any
- c now do inner points if there are any
- c
                 if(Maxdepth-k.gt.1) then
                   do k2=k+1,Maxdepth-1
- c
                     ttemp(k2) = ttemp(k2) +
       *              (mda(k2-1)*(tda(k2-1)-taa(k2-1))-
       *              mda(k2)*(tda(k2)-taa(k2)))
       *              *dt*recip_drF(k2)
- c
                    stemp(k2) = stemp(k2) +
       *              (mda(k2-1)*(sda(k2-1)-saa(k2-1))-
       *              mda(k2)*(sda(k2)-saa(k2)))
       *              *dt*recip_drF(k2)
- c
                   enddo
                 endif
                 ttemp(kmx+1) =  ttemp(kmx+1)+
       *                  (mda(kmx)*(tda(kmx)-taa(kmx)))*
       *                  dt*recip_drF(kmx+1)
- c
                 stemp(kmx+1) =  stemp(kmx+1)+
       *                  (mda(kmx)*(sda(kmx)-saa(kmx)))*
       *                  dt*recip_drF(kmx+1)
- c
- c set the environmental temp and salinity to equal new fields
+ C set the environmental temp and salinity to equal new fields
- c
                  do k2=1,maxdepth-1
                    taa(k2) = ttemp(k2+1)
                    saa(k2) = stemp(k2+1)
                  enddo
- c
- c end loop on number of time integration steps
+ C end loop on number of time integration steps
- c
               enddo
            ENDDO
       continue
- c
- c assume that it converged, so update the ta and sa with new fields
+ C assume that it converged, so update the ta and sa with new fields
- c
  c          if(i.gt.180.and.j.gt.200.and.i.lt.240) then
  c            write(*,*)"Converged ",i,j,k,maxdepth,ttemp(k+1),ta(i,k+1)
  c          endif
-Line 1126 
 c          endif
+Line 899 
 c          endif
  c          if(i.gt.180.and.j.gt.200.and.i.lt.240) then
  c            write(*,*)"convadj ",convadj(i,j,k2)
  c          endif
- c
- c see if nlopps messed up
+ C see if nlopps messed up
- c
              if(sa(i,k).gt.40..or.ta(i,k).lt.-4.0) then
                 problem = .true.
                 write(0,*)"t out of range at j after adjust",j
                 debug = .true.
              endif
- c
            enddo
- c
- c jump here if k = maxdepth or if level not unstable, go to next
+ C jump here if k = maxdepth or if level not unstable, go to next
- c profile point
+ C profile point
- c
      continue
- c
- c
+ C end loop on k, move on to next possible plume
- c end loop on k, move on to next possible plume
- c
        ENDDO
  continue
- c
- c i loop
+ C i loop
- c
        ENDDO
        END

 Legend:



Removed from v.1.10
 


changed lines


 
Added in v.1.11
 Legend:



Removed from v.1.10
 


changed lines


 
Added in v.1.11
-Removed from v.1.10
+Added in v.1.11

	ViewVC Help
Powered by ViewVC 1.1.22