C $Header: /u/gcmpack/MITgcm/pkg/opps/opps_calc.F,v 1.14 2016/03/10 20:57:11 jmc Exp $
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
C !INTERFACE: ======================================================
SUBROUTINE OPPS_CALC(
U tracerEnv,
I wVel,
I kMax, nTracer, nTracerInuse,
I I, J, bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *=====================================================================*
C | SUBROUTINE OPPS_CALC |
C | o Compute all OPPS fields defined in OPPS.h |
C *=====================================================================*
C | This subroutine is based on the routine 3dconvection.F |
C | by E. Skyllingstad (?) |
C | plenty of modifications to make it work: |
C | - removed many unused parameters and variables |
C | - turned everything (back) into 1D code |
C | - pass variables, that are originally in common blocks: |
C | maxDepth |
C | - pass vertical velocity, set in OPPS_INTERFACE |
C | - do not use convadj for now (whatever that is) |
C | - changed two .LT. 0 to .LE. 0 statements (because of possible |
C | division) |
C | - replaced statement function state1 by call to a real function |
C | - removed range check, actually moved it up to OPPS_INTERFACE |
C | - avoid division by zero: if (Wd.EQ.0) dt = ...1/Wd |
C | - cleaned-up debugging |
C | - replaced local dz and GridThickness by global drF |
C | - replaced 1/dz by 1*recip_drF |
C | - replaced 9.81 with gravity (=9.81) |
C | - added a lot of comments that relate code to equation in paper |
C | (Paluszkiewicz+Romea, 1997, Dynamics of Atmospheres and Oceans, |
C | 26, pp. 95-130) |
C | - included passive tracer support. This is the main change and may |
C | not improve the readability of the code because of the joint |
C | treatment of active (theta, salt) and passive tracers. The array |
C | tracerEnv(Nr,2+PTRACERS_num) contains |
C | theta = tracerEnv(:,1), |
C | salt = tracerEnv(:,2), and |
C | ptracers = tracerEnv(:,3:PTRACERS_num+2). |
C | All related array names have been changed accordingly, so that |
C | instead of Sd(Nr) and Td(Nr) (plume salinity and temperature), we |
C | have Pd(Nr,nTracer) (tracer in plume), with Sd(:) = Pd(:,2), |
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 \ev
IMPLICIT NONE
C !USES: ============================================================
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "OPPS.h"
#include "GRID.h"
C !INPUT PARAMETERS: ===================================================
C Routine arguments
C bi, bj :: array indices on which to apply calculations
C myTime :: Current time in simulation
INTEGER I, J, bi, bj, KMax, nTracer, nTracerInUse
INTEGER myThid, myIter
_RL myTime
_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 msgBuf :: Informational/error message buffer
CHARACTER*(MAX_LEN_MBUF) msgBuf
INTEGER K, K2, K2m1, K2p1, ktr
INTEGER ntime,nn,kmx,ic
INTEGER maxDepth
_RL wsqr,oldflux,newflux,entrainrate
_RL pmix
_RL D1, D2
_RL dz1,dz2
_RL radius,StartingFlux
_RL dtts,dt
C Arrays
_RL Paa(Nr,nTracer)
_RL wda(Nr), mda(Nr), pda(Nr,nTracer)
C
C Pd, Wd :: tracers, vertical velocity in plume
C Md :: plume mass flux (?)
C Ad :: fractional area covered by plume
C Dd :: density in plume
C De :: density of environment
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
C local initialization
C Copy some arrays
dtts = dTtracerLev(1)
C start k-loop
DO k=1,KMax-1
C initialize the plume T,S,density, and w velocity
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)
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)
CML compute vertical velocity at cell centers from GCM velocity
Wd(k)= - .5*(wVel(K)+wVel(K+1))
CML(
CML avoid division by zero
CML IF (Wd(K) .EQ. 0.D0) Wd(K) = 2.23e-16
CML)
C guess at initial top grid cell vertical velocity
CML Wd(k) = 0.03
C these estimates of initial plume velocity based on plume size and
C top grid cell water mass
c Wd(k) = 0.5*drF(k)/(dtts*FRACTIONAL_AREA)
c Wd(k) = 0.5*drF(k)/dtts
wsqr=Wd(k)*Wd(k)
PlumeEntrainment(k) = 0.0
#ifdef ALLOW_OPPS_DEBUG
IF ( OPPSdebugLevel.GE.debLevB ) THEN
WRITE(msgBuf,'(A,I3)')
& 'S/R OPPS_CALC: doing old lowerparcel', k
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
ENDIF
#endif /* ALLOW_OPPS_DEBUG */
radius=PlumeRadius
StartingFlux=radius*radius*Wd(k)*Dd(k)
oldflux=StartingFlux
dz2=DrF(k)
DO k2=k,KMax-1
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),
& i,j,k2+1,bi,bj,myThid)
De(k2+1)=D2
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 flux goes negative
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 find mass flux according to eq.(3) from paper by vertical integration
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
PlumeEntrainment(k2+1) = newflux/StartingFlux
IF(newflux.LE.0.0) then
#ifdef ALLOW_OPPS_DEBUG
IF ( OPPSdebugLevel.GE.debLevA ) THEN
WRITE(msgBuf,'(A,I3)')
& 'S/R OPPS_CALC: Plume entrained to zero at level ', k2
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
ENDIF
#endif /* ALLOW_OPPS_DEBUG */
maxdepth = k2
if(maxdepth.eq.k) goto 1000
goto 1
endif
C entrainment rate is basically a scaled mass flux dM/M
entrainrate = (newflux - oldflux)/newflux
oldflux = newflux
C mix var(s) are the average environmental values over the two grid levels
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 compute the density at this level for the buoyancy term in the
C vertical k.e. equation
Dd(k2+1)=STATE1(Pd(k2+1,2),Pd(k2+1,1),i,j,k2+1,bi,bj,myThid)
C next, solve for the vertical velocity k.e. using combined eq. (4)
C and eq (5) from the paper
#ifdef ALLOW_OPPS_DEBUG
IF ( OPPSdebugLevel.GE.debLevA ) THEN
WRITE(msgBuf,'(A,3E12.4,I3)')
& 'S/R OPPS_CALC: Dd,De,entr,k ',Dd(k2),De(k2),entrainrate,k2
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
ENDIF
#endif /* ALLOW_OPPS_DEBUG */
CML insert Eq. (4) into Eq. (5) to get something like this for wp^2
wsqr = wsqr - wsqr*abs(entrainrate)+ gravity*
& (dz1*(Dd(k2)-De(k2))/De(k2)
& +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1))
C if negative k.e. then plume has reached max depth, get out of loop
IF(wsqr.LE.0.0)then
maxdepth = k2
#ifdef ALLOW_OPPS_DEBUG
IF ( OPPSdebugLevel.GE.debLevA ) THEN
WRITE(msgBuf,'(A,I3)')
& 'S/R OPPS_CALC: Plume velocity went to zero at level ', k2
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
WRITE(msgBuf,'(A,4A14)')
& 'S/R OPPS_CALC: ', 'wsqr', 'entrainrate',
& '(Dd-De)/De up', '(Dd-De)/De do'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
WRITE(msgBuf,'(A,4E14.6)')
& 'S/R OPPS_CALC: ', wsqr, entrainrate,
& (Dd(k2)-De(k2))/De(k2), (Dd(k2+1)-De(k2+1))/De(k2+1)
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
ENDIF
#endif /* ALLOW_OPPS_DEBUG */
if(maxdepth.eq.k) goto 1000
goto 1
endif
Wd(k2+1)=sqrt(wsqr)
C compute a new radius based on the new mass flux at this grid level
C from Eq. (4)
radius=sqrt(newflux/(Wd(k2)*Dd(k2)))
ELSE
maxdepth=k2
if(maxdepth.eq.k) goto 1000
GOTO 1
ENDIF
ENDDO
C plume has reached the bottom
MaxDepth=kMax
1 CONTINUE
Ad(k)=FRACTIONAL_AREA
IC=0
C start iteration on fractional area, not used in OGCM implementation
DO IC=1,Max_ABE_Iterations
C next compute the mass flux beteen each grid box using the entrainment
Md(k)=Wd(k)*Ad(k)
DO k2=k+1,maxDepth
Md(k2)=Md(k)*PlumeEntrainment(k2)
#ifdef ALLOW_OPPS_DEBUG
IF ( OPPSdebugLevel.GE.debLevA ) THEN
WRITE(msgBuf,'(A,2E12.4,I3)')
& 'S/R OPPS_CALC: Md, Wd, and k are ',Md(k2),Wd(k2),k2
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
ENDIF
#endif /* ALLOW_OPPS_DEBUG */
ENDDO
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 center of grid cell, use weighted average to get boundary values
C use a timestep limited by the GCM model timestep and the maximum plume
C velocity (CFL criteria)
C calculate the weighted wd, td, and sd
dt = dtts
do k2=k,maxDepth-1
IF ( Wd(K2) .NE. 0. _d 0 ) dt = min(dt,drF(k2)/Wd(k2))
C time integration will be integer number of steps to get one
C gcm time step
ntime = nint(0.5*int(dtts/dt))
if(ntime.eq.0) then
ntime = 1
endif
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 cell. Eq. (20)
mda(k2) = (md(k2)*drF(k2)+md(k2+1)*drF(k2+1))/
* (drF(k2)+drF(k2+1))
wda(k2) = (wd(k2)*drF(k2)+wd(k2+1)*drF(k2+1))/
* (drF(k2)+drF(k2+1))
DO ktr = 1, nTracerInUse
Pda(k2,ktr) = Pd(k2,ktr)
Paa(k2,ktr) = tracerEnv(k2+1,ktr)
ENDDO
enddo
dt = min(dt,dtts)
#ifdef ALLOW_OPPS_DEBUG
IF ( OPPSdebugLevel.GE.debLevA ) THEN
WRITE(msgBuf,'(A,F14.4)')
& 'S/R OPPS_CALC: time step = ', dt
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
ENDIF
#endif /* ALLOW_OPPS_DEBUG */
DO ktr=1,nTracerInUse
Pda(maxdepth,ktr) = Pd(maxdepth,ktr)
ENDDO
kmx = maxdepth-1
do nn=1,ntime
C top point
DO ktr = 1,nTracerInUse
tracerEnv(k,ktr) = tracerEnv(k,ktr)-
& (mda(k)*(Pda(k,ktr)-Paa(k,ktr)))*dt*recip_drF(k)
ENDDO
C now do inner points if there are any
CML if(Maxdepth-k.gt.1) then
CML This if statement is superfluous
CML IF ( k .LT. Maxdepth-1 ) THEN
CML DO k2=k+1,Maxdepth-1
CML mda(maxDepth) = 0.
DO k2=k+1,kmx
k2m1 = max(k,k2-1)
k2p1 = max(k2+1,maxDepth)
DO ktr = 1,nTracerInUse
tracerEnv(k2,ktr) = tracerEnv(k2,ktr) +
& (mda(k2m1)*(Pda(k2m1,ktr)-Paa(k2m1,ktr))
& -mda(k2) *(Pda(k2,ktr) -Paa(k2,ktr)) )
& *dt*recip_drF(k2)
ENDDO
ENDDO
CML This if statement is superfluous
CML ENDIF
C bottom point
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 set the environmental temp and salinity to equal new fields
DO ktr=1,nTracerInUse
DO k2=1,kmx
paa(k2,ktr) = tracerEnv(k2+1,ktr)
ENDDO
ENDDO
C end loop on number of time integration steps
enddo
ENDDO
C count convection event in this grid cell
OPPSconvectCount(I,J,K,bi,bj) =
& OPPSconvectCount(I,J,K,bi,bj) + 1. _d 0
C jump here if k = maxdepth or if level not unstable, go to next
C profile point
1000 continue
C end of k-loop
ENDDO
C-- End IF (DIFFERENT_MULTIPLE)
ENDIF
RETURN
END
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
C | Calculates rho(S,T,p)
C | It is absolutely necessary to compute
C | the full rho and not sigma=rho-rhoConst, because
C | density is used as a scale factor for fluxes and velocities
C *===============================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_NONHYDROSTATIC
# include "NH_VARS.h"
#endif /* ALLOW_NONHYDROSTATIC */
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
INTEGER i, j, kRef, bi, bj, myThid
_RL tLoc, sLoc
C !LOCAL VARIABLES:
C == Local variables ==
_RL rhoLoc
_RL pLoc
CMLC estimate pressure from depth at cell centers
CML mtoSI = gravity*rhoConst
CML pLoc = ABS(rC(kRef))*mtoSI
IF ( usingZCoords ) THEN
C in Z coordinates the pressure is rho0 * (hydrostatic) Potential
#ifdef ALLOW_NONHYDROSTATIC
IF ( selectP_inEOS_Zc.EQ.3 ) THEN
pLoc = rhoConst*( totPhiHyd(i,j,kRef,bi,bj)
& + phi_nh(i,j,kRef,bi,bj)
& + phiRef(2*kRef)
& )*maskC(i,j,kRef,bi,bj)
ELSEIF ( selectP_inEOS_Zc.EQ.2 ) THEN
#else /* ALLOW_NONHYDROSTATIC */
IF ( selectP_inEOS_Zc.EQ.2 ) THEN
#endif /* ALLOW_NONHYDROSTATIC */
C----------
C NOTE: For now, totPhiHyd only contains the Potential anomaly
C since PhiRef has not (yet) been added in S/R DIAGS_PHI_HYD
C----------
pLoc = rhoConst*( totPhiHyd(i,j,kRef,bi,bj)
& + phiRef(2*kRef)
& )*maskC(i,j,kRef,bi,bj)
c ELSEIF ( selectP_inEOS_Zc.EQ.1 ) THEN
C note: for the case selectP_inEOS_Zc=0, also use pRef4EOS (set to
C rhoConst*phiRef(2*k) ) to reproduce same previous machine truncation
ELSEIF ( selectP_inEOS_Zc.LE.1 ) THEN
pLoc = pRef4EOS(kRef)*maskC(i,j,kRef,bi,bj)
ELSE
pLoc = rhoConst*phiRef(2*kRef)*maskC(i,j,kRef,bi,bj)
ENDIF
ELSEIF ( usingPCoords ) THEN
C in P coordinates the pressure is just the coordinate of the tracer point
pLoc = rC(kRef)* maskC(i,j,kRef,bi,bj)
ENDIF
CALL FIND_RHO_SCALAR( tLoc, sLoc, pLoc, rhoLoc, myThid )
STATE1 = rhoLoc
#endif /* ALLOW_OPPS */
RETURN
END
#ifdef OPPS_ORGCODE
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
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 In the POP code. The subroutine is called nlopps (long story to explain why).
C I've attached the version of lopps that we've been using in the simulations.
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 does need "kmp" which is why the boundc common is included. For massively
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)
parameter (imt = 361 , jmt = 301 , km = 30 )
C Nlopps: E. Skyllingstad and T. Paluszkiewicz
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"
c include "plume_fast.h"
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
INTEGER i,j,k
Clfh
integer is,ie,k2
Clfh
REAL D1,D2,state1,Density
REAL dz1,dz2
REAL radius,StartingFlux
real ttemp(km),stemp(km),taa(km),saa(km)
real wda(km),tda(km),sda(km),mda(km)
real dtts,dt,sumo,sumn
integer ntime,nn,kmx,ic
LOGICAL debug,done
INTEGER MAX_ABE_ITERATIONS
PARAMETER(MAX_ABE_ITERATIONS=1)
REAL PlumeRadius
REAL STABILITY_THRESHOLD
REAL FRACTIONAL_AREA
REAL MAX_FRACTIONAL_AREA
REAL VERTICAL_VELOCITY
REAL ENTRAINMENT_RATE
REAL e2
PARAMETER ( PlumeRadius = 100.D0 )
PARAMETER ( STABILITY_THRESHOLD = -1.E-4 )
PARAMETER ( FRACTIONAL_AREA = .1E0 )
PARAMETER ( MAX_FRACTIONAL_AREA = .8E0 )
PARAMETER ( VERTICAL_VELOCITY = .02E0 )
PARAMETER ( ENTRAINMENT_RATE = -.05E0 )
PARAMETER ( e2 = 2.E0*ENTRAINMENT_RATE )
! Arrays.
REAL Ad(km),Sd(km),Td(km),Wd(km),Dd(km),Md(km)
REAL Se(km),Te(km),We(km),De(km)
REAL PlumeEntrainment(km)
REAL GridThickness(km)
C input kmp through a common block
common / boundc / wsx(imt,jmt),wsy(imt,jmt),hfs(imt,jmt),
1 ple(imt,jmt),kmp(imt,jmt),kmq(imt,jmt)
cwmseas
& ,wsx1(imt,jmt),wsy1(imt,jmt)
1 ,wsx2(imt,jmt),wsy2(imt,jmt)
C input the variables through a common
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 otherwise it is repetive
C griddz is initialize by call to setupgrid
dtts = 2400
do k=1,km
dz(k) = 0.01*gcmdz(k)
enddo
do k=1,km
GridThickness(k) = dz(k)
enddo
C modified to loop over slab
DO i=is,ie
numgridpoints=kmp(i,j)
C go to next column if only 1 grid point or on land
if(numgridpoints.le.1) goto 1100
C loop over depth
debug = .false.
C first save copy of initial profile
DO k=1,NumGridPoints
stemp(k)=sa(i,k)
ttemp(k)=ta(i,k)
C do a check of t and s range, if out of bounds set flag
if(problem) then
write(0,*)"Code in trouble before this nlopps call"
return
endif
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
debug = .true.
return
endif
ENDDO
if(debug) then
write(*,*)"T and S Profile at ",i,j
write(*,*)(ta(i,k),sa(i,k),k=1,NumGridPoints)
endif
DO k=1,NumGridPoints-1
C initialize the plume T,S,density, and w velocity
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 guess at initial top grid cell vertical velocity
Wd(k) = 0.03
C these estimates of initial plume velocity based on plume size and
C top grid cell water mass
c Wd(k) = 0.5*dz(k)/(dtts*FRACTIONAL_AREA)
c Wd(k) = 0.5*dz(k)/dtts
wsqr=Wd(k)*Wd(k)
PlumeEntrainment(k) = 0.0
if(debug) write(0,*) 'Doing old lowerparcel'
radius=PlumeRadius
StartingFlux=radius*radius*Wd(k)*Dd(k)
oldflux=StartingFlux
dz2=GridThickness(k)
DO k2=k,NumGridPoints-1
D1=state1(Sd(k2),Td(k2),k2+1)
D2=state1(stemp(k2+1),ttemp(k2+1),k2+1)
De(k2+1)=D2
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 flux goes negative
IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN
dz1=dz2
dz2=GridThickness(k2+1)
C define mass flux according to eq. 4 from paper
newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*0.50*
. (dz1+dz2)
PlumeEntrainment(k2+1) = newflux/StartingFlux
IF(newflux.LT.0.0) then
if(debug) then
write(0,*)"Plume entrained to zero at ",k2
endif
maxdepth = k2
if(maxdepth.eq.k) goto 1000
goto 1
endif
C entrainment rate is basically a scaled mass flux dM/M
entrainrate = (newflux - oldflux)/newflux
oldflux = newflux
C mix var(s) are the average environmental values over the two grid levels
smix=(dz1*stemp(k2)+dz2*stemp(k2+1))/(dz1+dz2)
thmix=(dz1*ttemp(k2)+dz2*ttemp(k2+1))/(dz1+dz2)
C first compute the new salinity and temperature for this level
C using equations 3.6 and 3.7 from the paper
sd(k2+1)=sd(k2) - entrainrate*(smix - sd(k2))
td(k2+1)=td(k2) - entrainrate*(thmix - td(k2))
C compute the density at this level for the buoyancy term in the
C vertical k.e. equation
Dd(k2+1)=state1(Sd(k2+1),Td(k2+1),k2+1)
C next, solve for the vertical velocity k.e. using combined eq. 4
C and eq 5 from the paper
if(debug) then
write(0,*)"Dd,De,entr,k ",Dd(k2),De(k2),entrainrate,k2
endif
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 if negative k.e. then plume has reached max depth, get out of loop
IF(wsqr.LT.0.0)then
maxdepth = k2
if(debug) then
write(0,*)"Plume velocity went to zero at ",k2
endif
if(maxdepth.eq.k) goto 1000
goto 1
endif
Wd(k2+1)=sqrt(wsqr)
C compute a new radius based on the new mass flux at this grid level
radius=sqrt(newflux/(Wd(k2)*Dd(k2)))
ELSE
maxdepth=k2
if(maxdepth.eq.k) goto 1000
GOTO 1
ENDIF
ENDDO
C plume has reached the bottom
MaxDepth=NumGridPoints
1 continue
Ad(k)=FRACTIONAL_AREA
IC=0
C start iteration on fractional area, not used in OGCM implementation
DO IC=1,Max_ABE_Iterations
C next compute the mass flux beteen each grid box using the entrainment
92 continue
Md(k)=Wd(k)*Ad(k)
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 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 center of grid cell, use weighted average to get boundary values
C use a timestep limited by the GCM model timestep and the maximum plume
C velocity (CFL criteria)
C calculate the weighted wd, td, and sd
dt = dtts
do k2=k,maxdepth-1
dt = min(dt,dz(k2)/wd(k2))
C time integration will be integer number of steps to get one
C gcm time step
ntime = nint(0.5*int(dtts/dt))
if(ntime.eq.0) then
ntime = 1
endif
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 cell.
mda(k2) = (md(k2)*dz(k2)+md(k2+1)*dz(k2+1))/
* (dz(k2)+dz(k2+1))
wda(k2) = (wd(k2)*dz(k2)+wd(k2+1)*dz(k2+1))/
* (dz(k2)+dz(k2+1))
tda(k2) = td(k2)
sda(k2) = sd(k2)
taa(k2) = ttemp(k2+1)
saa(k2) = stemp(k2+1)
enddo
dt = min(dt,dtts)
if(debug) then
write(0,*)"Time step is ", dt
endif
tda(maxdepth) = td(maxdepth)
sda(maxdepth) = sd(maxdepth)
C do top and bottom points first
kmx = maxdepth-1
do nn=1,ntime
ttemp(k) = ttemp(k)-
* (mda(k)*(tda(k)-taa(k)))*dt*recip_drF(k)
stemp(k) = stemp(k)-
* (mda(k)*(sda(k)-saa(k)))*dt*recip_drF(k)
C now do inner points if there are any
if(Maxdepth-k.gt.1) then
do k2=k+1,Maxdepth-1
ttemp(k2) = ttemp(k2) +
* (mda(k2-1)*(tda(k2-1)-taa(k2-1))-
* mda(k2)*(tda(k2)-taa(k2)))
* *dt*recip_drF(k2)
stemp(k2) = stemp(k2) +
* (mda(k2-1)*(sda(k2-1)-saa(k2-1))-
* mda(k2)*(sda(k2)-saa(k2)))
* *dt*recip_drF(k2)
enddo
endif
ttemp(kmx+1) = ttemp(kmx+1)+
* (mda(kmx)*(tda(kmx)-taa(kmx)))*
* dt*recip_drF(kmx+1)
stemp(kmx+1) = stemp(kmx+1)+
* (mda(kmx)*(sda(kmx)-saa(kmx)))*
* dt*recip_drF(kmx+1)
C set the environmental temp and salinity to equal new fields
do k2=1,maxdepth-1
taa(k2) = ttemp(k2+1)
saa(k2) = stemp(k2+1)
enddo
C end loop on number of time integration steps
enddo
ENDDO
999 continue
C assume that it converged, so update the ta and sa with new fields
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
do k2=k,maxdepth
convadj(i,j,k2) = convadj(i,j,k2) + (ttemp(k2)-
* ta(i,k2))
sa(i,k2) = stemp(k2)
ta(i,k2) = ttemp(k2)
c if(i.gt.180.and.j.gt.200.and.i.lt.240) then
c write(*,*)"convadj ",convadj(i,j,k2)
c endif
C see if nlopps messed up
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
enddo
C jump here if k = maxdepth or if level not unstable, go to next
C profile point
1000 continue
C end loop on k, move on to next possible plume
ENDDO
1100 continue
C i loop
ENDDO
END
#endif /* OPPS_ORGCODE */