pkg/opps/opps_calc.F

893df04db0 Mart*0001 #include "OPPS_OPTIONS.h"
6580feb7eb Jean*0002 #undef OPPS_ORGCODE
                0003
                0004 C--  File opps_calc.F:
                0005 C--   Contents
                0006 C--   o OPPS_CALC
                0007 C--   o STATE1
                0008 C--   o NLOPPS
893df04db0 Mart*0009
6580feb7eb Jean*0010 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
893df04db0 Mart*0011 CBOP
                0012 C !ROUTINE: OPPS_CALC
                0013
                0014 C !INTERFACE: ======================================================
a6de755a30 Jean*0015       SUBROUTINE OPPS_CALC(
893df04db0 Mart*0016      U     tracerEnv,
016b84c482 Mart*0017      O     OPPSconvectCount,
3daafce20b Jean*0018      I     wVel,
                0019      I     kMax, nTracer, nTracerInuse,
893df04db0 Mart*0020      I     I, J, bi, bj, myTime, myIter, myThid )
                0021
                0022 C !DESCRIPTION: \bv
6580feb7eb Jean*0023 C     *=====================================================================*
893df04db0 Mart*0024 C     | SUBROUTINE OPPS_CALC                                                |
                0025 C     | o Compute all OPPS fields defined in OPPS.h                         |
6580feb7eb Jean*0026 C     *=====================================================================*
893df04db0 Mart*0027 C     | This subroutine is based on the routine 3dconvection.F              |
                0028 C     | by E. Skyllingstad (?)                                              |
                0029 C     | plenty of modifications to make it work:                            |
                0030 C     | - removed many unused parameters and variables                      |
                0031 C     | - turned everything (back) into 1D code                             |
29188f9691 Jean*0032 C     | - pass variables, that are originally in common blocks:             |
893df04db0 Mart*0033 C     |   maxDepth                                                          |
                0034 C     | - pass vertical velocity, set in OPPS_INTERFACE                     |
                0035 C     | - do not use convadj for now (whatever that is)                     |
                0036 C     | - changed two .LT. 0 to .LE. 0 statements (because of possible      |
                0037 C     |   division)                                                         |
                0038 C     | - replaced statement function state1 by call to a real function     |
                0039 C     | - removed range check, actually moved it up to OPPS_INTERFACE       |
                0040 C     | - avoid division by zero: if (Wd.EQ.0) dt = ...1/Wd                 |
                0041 C     | - cleaned-up debugging                                              |
                0042 C     | - replaced local dz and GridThickness by global drF                 |
                0043 C     | - replaced 1/dz by 1*recip_drF                                      |
                0044 C     | - replaced 9.81 with gravity (=9.81)                                |
                0045 C     | - added a lot of comments that relate code to equation in paper     |
                0046 C     |   (Paluszkiewicz+Romea, 1997, Dynamics of Atmospheres and Oceans,   |
                0047 C     |   26, pp. 95-130)                                                   |
                0048 C     | - included passive tracer support. This is the main change and may  |
                0049 C     |   not improve the readability of the code because of the joint      |
                0050 C     |   treatment of active (theta, salt) and passive tracers. The array  |
                0051 C     |   tracerEnv(Nr,2+PTRACERS_num) contains                             |
                0052 C     |   theta    = tracerEnv(:,1),                                        |
3daafce20b Jean*0053 C     |   salt     = tracerEnv(:,2), and                                    |
893df04db0 Mart*0054 C     |   ptracers = tracerEnv(:,3:PTRACERS_num+2).                         |
                0055 C     |   All related array names have been changed accordingly, so that    |
                0056 C     |   instead of Sd(Nr) and Td(Nr) (plume salinity and temperature), we |
                0057 C     |   have Pd(Nr,nTracer) (tracer in plume), with Sd(:) = Pd(:,2),      |
                0058 C     |   Td(:) = Pd(:,1), etc.                                             |
                0059 C     | o TODO:                                                             |
                0060 C     |   clean up the logic of the vertical loops and get rid off the      |
                0061 C     |   GOTO statements                                                   |
6580feb7eb Jean*0062 C     *=====================================================================*
893df04db0 Mart*0063 C \ev
                0064
6580feb7eb Jean*0065       IMPLICIT NONE
                0066
016b84c482 Mart*0067 C !USES: ==============================================================
893df04db0 Mart*0068 #include "SIZE.h"
                0069 #include "EEPARAMS.h"
                0070 #include "PARAMS.h"
                0071 #include "OPPS.h"
                0072 #include "GRID.h"
                0073
016b84c482 Mart*0074 C !INPUT/OUTPUT PARAMETERS: ============================================
6580feb7eb Jean*0075 C Routine arguments
016b84c482 Mart*0076 C     OPPSconvectCount :: counter for freqency of convection events
6580feb7eb Jean*0077 C     bi, bj :: array indices on which to apply calculations
                0078 C     myTime :: Current time in simulation
016b84c482 Mart*0079       INTEGER KMax, nTracer, nTracerInUse
                0080       _RL tracerEnv(Nr,nTracer)
                0081       _RL OPPSconvectCount(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0082       _RL wVel(Nr)
                0083       INTEGER I, J, bi, bj
893df04db0 Mart*0084       _RL     myTime
016b84c482 Mart*0085       INTEGER myThid, myIter
893df04db0 Mart*0086
                0087 #ifdef ALLOW_OPPS
6580feb7eb Jean*0088 C !FUNCTIONS: ==========================================================
                0089 c     EXTERNAL DIFFERENT_MULTIPLE
                0090 c     LOGICAL  DIFFERENT_MULTIPLE
                0091       _RL STATE1
                0092
893df04db0 Mart*0093 C !LOCAL VARIABLES: ====================================================
6580feb7eb Jean*0094 C Local constants
                0095 C     msgBuf      :: Informational/error message buffer
893df04db0 Mart*0096       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0097       INTEGER K, K2, K2m1, K2p1, ktr
                0098       INTEGER ntime,nn,kmx,ic
                0099       INTEGER maxDepth
                0100
                0101       _RL wsqr,oldflux,newflux,entrainrate
                0102       _RL pmix
6580feb7eb Jean*0103       _RL D1, D2
893df04db0 Mart*0104       _RL dz1,dz2
                0105       _RL radius,StartingFlux
                0106       _RL dtts,dt
                0107 C     Arrays
                0108       _RL Paa(Nr,nTracer)
                0109       _RL wda(Nr), mda(Nr), pda(Nr,nTracer)
                0110 C
6580feb7eb Jean*0111 C     Pd, Wd           :: tracers, vertical velocity in plume
                0112 C     Md               :: plume mass flux (?)
                0113 C     Ad               :: fractional area covered by plume
                0114 C     Dd               :: density in plume
                0115 C     De               :: density of environment
                0116 C     PlumeEntrainment ::
893df04db0 Mart*0117       _RL Ad(Nr),Wd(Nr),Dd(Nr),Md(Nr)
                0118       _RL De(Nr)
                0119       _RL PlumeEntrainment(Nr)
                0120       _RL Pd(Nr,nTracer)
                0121 CEOP
                0122
                0123 C--   Check to see if should convect now
94a46dfe0d Jean*0124 C      IF ( DIFFERENT_MULTIPLE(cAdjFreq,myTime,deltaTClock) ) THEN
893df04db0 Mart*0125       IF ( .true. ) THEN
                0126 C     local initialization
                0127
                0128 C     Copy some arrays
dfc17c9c63 Jean*0129       dtts = dTtracerLev(1)
6580feb7eb Jean*0130
893df04db0 Mart*0131 C     start k-loop
                0132
                0133       DO k=1,KMax-1
6580feb7eb Jean*0134
                0135 C initialize the plume T,S,density, and w velocity
                0136
893df04db0 Mart*0137        DO ktr=1,nTracerInUse
                0138         Pd(k,ktr) = tracerEnv(k,ktr)
                0139        ENDDO
6580feb7eb Jean*0140        Dd(k)=STATE1(Pd(k,2),Pd(k,1),i,j,k,bi,bj,myThid)
893df04db0 Mart*0141        De(k)=Dd(k)
                0142 CML       print *, 'ml-opps:', i,j,k,tracerEnv(k,2),tracerEnv(k,1),
                0143 CML     &      Dd(k),Pd(k,1),Pd(k,2)
                0144 CML compute vertical velocity at cell centers from GCM velocity
                0145        Wd(k)= - .5*(wVel(K)+wVel(K+1))
                0146 CML(
98bf704dd5 Jean*0147 CML    avoid division by zero
893df04db0 Mart*0148 CML       IF (Wd(K) .EQ. 0.D0) Wd(K) = 2.23e-16
                0149 CML)
6580feb7eb Jean*0150
                0151 C guess at initial top grid cell vertical velocity
                0152
893df04db0 Mart*0153 CML          Wd(k) = 0.03
6580feb7eb Jean*0154
                0155 C these estimates of initial plume velocity based on plume size and
                0156 C top grid cell water mass
                0157
893df04db0 Mart*0158 c          Wd(k) = 0.5*drF(k)/(dtts*FRACTIONAL_AREA)
                0159 c          Wd(k) = 0.5*drF(k)/dtts
6580feb7eb Jean*0160
893df04db0 Mart*0161        wsqr=Wd(k)*Wd(k)
                0162        PlumeEntrainment(k) = 0.0
6580feb7eb Jean*0163
893df04db0 Mart*0164 #ifdef ALLOW_OPPS_DEBUG
                0165        IF ( OPPSdebugLevel.GE.debLevB ) THEN
                0166         WRITE(msgBuf,'(A,I3)')
3daafce20b Jean*0167      &       'S/R OPPS_CALC: doing old lowerparcel', k
893df04db0 Mart*0168         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0169      &       SQUEEZE_RIGHT , 1)
                0170        ENDIF
                0171 #endif /* ALLOW_OPPS_DEBUG */
                0172        radius=PlumeRadius
                0173        StartingFlux=radius*radius*Wd(k)*Dd(k)
                0174        oldflux=StartingFlux
3daafce20b Jean*0175
893df04db0 Mart*0176        dz2=DrF(k)
                0177        DO k2=k,KMax-1
6580feb7eb Jean*0178         D1=STATE1( Pd(k2,2), Pd(k2,1),i,j,k2+1,bi,bj,myThid)
                0179         D2=STATE1( tracerEnv(k2+1,2), tracerEnv(k2+1,1),
893df04db0 Mart*0180      &                                i,j,k2+1,bi,bj,myThid)
                0181         De(k2+1)=D2
6580feb7eb Jean*0182
                0183 C To start downward, parcel has to initially be heavier than environment
                0184 C but after it has started moving, we continue plume until plume tke or
                0185 C flux goes negative
                0186
893df04db0 Mart*0187 CML     &       _hFacC(i,j,k-1,bi,bj)
                0188 CML     &       *_hFacC(i,j,k,bi,bj) .GT. 0.
                0189 CML     &  .AND.
                0190         IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN
                0191          dz1=dz2
                0192          dz2=DrF(k2+1)
6580feb7eb Jean*0193
893df04db0 Mart*0194 C     find mass flux according to eq.(3) from paper by vertical integration
6580feb7eb Jean*0195
893df04db0 Mart*0196          newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*
                0197      &        .5*(dz1+dz2)
                0198 CML         print *, 'ml-opps:', i,j,k,oldflux,newflux,e2,radius,
                0199 CML     &        Wd(k2),Dd(k2),Pd(k2,1),Pd(k2,2),dz1,dz2
6580feb7eb Jean*0200
893df04db0 Mart*0201          PlumeEntrainment(k2+1) = newflux/StartingFlux
6580feb7eb Jean*0202
893df04db0 Mart*0203          IF(newflux.LE.0.0) then
                0204 #ifdef ALLOW_OPPS_DEBUG
                0205           IF ( OPPSdebugLevel.GE.debLevA ) THEN
3daafce20b Jean*0206            WRITE(msgBuf,'(A,I3)')
893df04db0 Mart*0207      &          'S/R OPPS_CALC: Plume entrained to zero at level ', k2
                0208            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0209      &          SQUEEZE_RIGHT , 1)
                0210           ENDIF
                0211 #endif /* ALLOW_OPPS_DEBUG */
                0212           maxdepth = k2
                0213           if(maxdepth.eq.k) goto 1000
                0214           goto 1
                0215          endif
6580feb7eb Jean*0216
                0217 C entrainment rate is basically a scaled mass flux dM/M
                0218
893df04db0 Mart*0219          entrainrate = (newflux - oldflux)/newflux
                0220          oldflux = newflux
6580feb7eb Jean*0221
                0222 C mix var(s) are the average environmental values over the two grid levels
                0223
893df04db0 Mart*0224          DO ktr=1,nTracerInUse
                0225           pmix=(dz1*tracerEnv(k2,ktr)+dz2*tracerEnv(k2+1,ktr))
                0226      &         /(dz1+dz2)
3daafce20b Jean*0227           Pd(k2+1,ktr)=Pd(k2,ktr)
893df04db0 Mart*0228      &         - entrainrate*(pmix - Pd(k2,ktr))
                0229          ENDDO
6580feb7eb Jean*0230
                0231 C compute the density at this level for the buoyancy term in the
                0232 C vertical k.e. equation
                0233
                0234          Dd(k2+1)=STATE1(Pd(k2+1,2),Pd(k2+1,1),i,j,k2+1,bi,bj,myThid)
                0235
                0236 C next, solve for the vertical velocity k.e. using combined eq. (4)
                0237 C and eq (5) from the paper
                0238
893df04db0 Mart*0239 #ifdef ALLOW_OPPS_DEBUG
                0240          IF ( OPPSdebugLevel.GE.debLevA ) THEN
3daafce20b Jean*0241           WRITE(msgBuf,'(A,3E12.4,I3)')
893df04db0 Mart*0242      &    'S/R OPPS_CALC: Dd,De,entr,k ',Dd(k2),De(k2),entrainrate,k2
                0243           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0244      &         SQUEEZE_RIGHT , 1)
                0245          ENDIF
                0246 #endif /* ALLOW_OPPS_DEBUG */
                0247 CML   insert Eq. (4) into Eq. (5) to get something like this for wp^2
                0248          wsqr = wsqr - wsqr*abs(entrainrate)+ gravity*
                0249      &        (dz1*(Dd(k2)-De(k2))/De(k2)
                0250      &        +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1))
6580feb7eb Jean*0251
                0252 C if negative k.e. then plume has reached max depth, get out of loop
                0253
893df04db0 Mart*0254          IF(wsqr.LE.0.0)then
                0255           maxdepth = k2
                0256 #ifdef ALLOW_OPPS_DEBUG
                0257           IF ( OPPSdebugLevel.GE.debLevA ) THEN
3daafce20b Jean*0258            WRITE(msgBuf,'(A,I3)')
893df04db0 Mart*0259      &     'S/R OPPS_CALC: Plume velocity went to zero at level ', k2
                0260            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0261      &          SQUEEZE_RIGHT , 1)
3daafce20b Jean*0262            WRITE(msgBuf,'(A,4A14)')
                0263      &          'S/R OPPS_CALC: ', 'wsqr', 'entrainrate',
893df04db0 Mart*0264      &          '(Dd-De)/De up', '(Dd-De)/De do'
                0265            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0266      &          SQUEEZE_RIGHT , 1)
3daafce20b Jean*0267            WRITE(msgBuf,'(A,4E14.6)')
                0268      &          'S/R OPPS_CALC: ', wsqr, entrainrate,
893df04db0 Mart*0269      &          (Dd(k2)-De(k2))/De(k2), (Dd(k2+1)-De(k2+1))/De(k2+1)
                0270            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0271      &          SQUEEZE_RIGHT , 1)
                0272           ENDIF
                0273 #endif /* ALLOW_OPPS_DEBUG */
                0274           if(maxdepth.eq.k) goto 1000
                0275           goto 1
                0276          endif
                0277          Wd(k2+1)=sqrt(wsqr)
6580feb7eb Jean*0278
893df04db0 Mart*0279 C     compute a new radius based on the new mass flux at this grid level
                0280 C     from Eq. (4)
6580feb7eb Jean*0281
893df04db0 Mart*0282          radius=sqrt(newflux/(Wd(k2)*Dd(k2)))
                0283         ELSE
                0284          maxdepth=k2
                0285          if(maxdepth.eq.k) goto 1000
                0286          GOTO 1
                0287         ENDIF
                0288        ENDDO
6580feb7eb Jean*0289
                0290 C plume has reached the bottom
                0291
893df04db0 Mart*0292        MaxDepth=kMax
6580feb7eb Jean*0293
893df04db0 Mart*0294  1     CONTINUE
6580feb7eb Jean*0295
893df04db0 Mart*0296        Ad(k)=FRACTIONAL_AREA
                0297        IC=0
6580feb7eb Jean*0298
                0299 C start iteration on fractional area, not used in OGCM implementation
                0300
893df04db0 Mart*0301        DO IC=1,Max_ABE_Iterations
6580feb7eb Jean*0302
                0303 C next compute the mass flux beteen each grid box using the entrainment
                0304
893df04db0 Mart*0305         Md(k)=Wd(k)*Ad(k)
6580feb7eb Jean*0306
893df04db0 Mart*0307         DO k2=k+1,maxDepth
                0308          Md(k2)=Md(k)*PlumeEntrainment(k2)
                0309 #ifdef ALLOW_OPPS_DEBUG
                0310          IF ( OPPSdebugLevel.GE.debLevA ) THEN
3daafce20b Jean*0311           WRITE(msgBuf,'(A,2E12.4,I3)')
893df04db0 Mart*0312      &         'S/R OPPS_CALC: Md, Wd, and  k are ',Md(k2),Wd(k2),k2
                0313           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0314      &         SQUEEZE_RIGHT , 1)
                0315          ENDIF
                0316 #endif /* ALLOW_OPPS_DEBUG */
                0317         ENDDO
6580feb7eb Jean*0318
                0319 C Now move on to calculate new temperature using flux from
                0320 C Td, Sd, Wd, ta, sa, and we. Values for these variables are at
                0321 C center of grid cell, use weighted average to get boundary values
                0322
                0323 C use a timestep limited by the GCM model timestep and the maximum plume
                0324 C velocity (CFL criteria)
                0325
                0326 C calculate the weighted wd, td, and sd
893df04db0 Mart*0327         dt = dtts
                0328         do k2=k,maxDepth-1
                0329          IF ( Wd(K2) .NE. 0. _d 0 ) dt = min(dt,drF(k2)/Wd(k2))
6580feb7eb Jean*0330
                0331 C time integration will be integer number of steps to get one
                0332 C gcm time step
                0333
893df04db0 Mart*0334          ntime = nint(0.5*int(dtts/dt))
                0335          if(ntime.eq.0) then
                0336           ntime = 1
                0337          endif
6580feb7eb Jean*0338
                0339 C make sure area weighted vertical velocities match; in other words
                0340 C make sure mass in equals mass out at the intersection of each grid
                0341 C cell. Eq. (20)
                0342
893df04db0 Mart*0343          mda(k2) = (md(k2)*drF(k2)+md(k2+1)*drF(k2+1))/
                0344      *        (drF(k2)+drF(k2+1))
6580feb7eb Jean*0345
893df04db0 Mart*0346          wda(k2) = (wd(k2)*drF(k2)+wd(k2+1)*drF(k2+1))/
                0347      *        (drF(k2)+drF(k2+1))
6580feb7eb Jean*0348
893df04db0 Mart*0349          DO ktr = 1, nTracerInUse
                0350           Pda(k2,ktr) = Pd(k2,ktr)
                0351           Paa(k2,ktr) = tracerEnv(k2+1,ktr)
98bf704dd5 Jean*0352          ENDDO
6580feb7eb Jean*0353
893df04db0 Mart*0354         enddo
                0355         dt = min(dt,dtts)
                0356 #ifdef ALLOW_OPPS_DEBUG
                0357         IF ( OPPSdebugLevel.GE.debLevA ) THEN
3daafce20b Jean*0358          WRITE(msgBuf,'(A,F14.4)')
893df04db0 Mart*0359      &        'S/R OPPS_CALC: time step = ', dt
                0360          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0361      &        SQUEEZE_RIGHT , 1)
                0362         ENDIF
                0363 #endif /* ALLOW_OPPS_DEBUG */
                0364         DO ktr=1,nTracerInUse
                0365          Pda(maxdepth,ktr) = Pd(maxdepth,ktr)
98bf704dd5 Jean*0366         ENDDO
6580feb7eb Jean*0367
893df04db0 Mart*0368         kmx = maxdepth-1
                0369         do nn=1,ntime
6580feb7eb Jean*0370
3daafce20b Jean*0371 C     top point
6580feb7eb Jean*0372
893df04db0 Mart*0373          DO ktr = 1,nTracerInUse
                0374           tracerEnv(k,ktr) =  tracerEnv(k,ktr)-
                0375      &        (mda(k)*(Pda(k,ktr)-Paa(k,ktr)))*dt*recip_drF(k)
98bf704dd5 Jean*0376          ENDDO
6580feb7eb Jean*0377
                0378 C now do inner points if there are any
                0379
893df04db0 Mart*0380 CML         if(Maxdepth-k.gt.1) then
3daafce20b Jean*0381 CML    This if statement is superfluous
893df04db0 Mart*0382 CML         IF ( k .LT. Maxdepth-1 ) THEN
                0383 CML         DO k2=k+1,Maxdepth-1
                0384 CML         mda(maxDepth) = 0.
                0385          DO k2=k+1,kmx
                0386           k2m1 = max(k,k2-1)
                0387           k2p1 = max(k2+1,maxDepth)
6580feb7eb Jean*0388
893df04db0 Mart*0389            DO ktr = 1,nTracerInUse
                0390             tracerEnv(k2,ktr) = tracerEnv(k2,ktr) +
                0391      &           (mda(k2m1)*(Pda(k2m1,ktr)-Paa(k2m1,ktr))
                0392      &           -mda(k2)  *(Pda(k2,ktr)  -Paa(k2,ktr))  )
                0393      &           *dt*recip_drF(k2)
3daafce20b Jean*0394            ENDDO
893df04db0 Mart*0395           ENDDO
3daafce20b Jean*0396 CML    This if statement is superfluous
893df04db0 Mart*0397 CML         ENDIF
6580feb7eb Jean*0398
893df04db0 Mart*0399 C     bottom point
6580feb7eb Jean*0400
893df04db0 Mart*0401          DO ktr=1,nTracerInUse
                0402           tracerEnv(kmx+1,ktr) =  tracerEnv(kmx+1,ktr)+
                0403      &        mda(kmx)*(Pda(kmx,ktr)-Paa(kmx,ktr))*dt*recip_drF(kmx+1)
                0404          ENDDO
6580feb7eb Jean*0405
                0406 C     set the environmental temp and salinity to equal new fields
                0407
893df04db0 Mart*0408          DO ktr=1,nTracerInUse
                0409           DO k2=1,kmx
                0410            paa(k2,ktr) = tracerEnv(k2+1,ktr)
                0411           ENDDO
98bf704dd5 Jean*0412          ENDDO
6580feb7eb Jean*0413
                0414 C end loop on number of time integration steps
                0415
893df04db0 Mart*0416         enddo
                0417        ENDDO
6580feb7eb Jean*0418
893df04db0 Mart*0419 C     count convection event in this grid cell
6580feb7eb Jean*0420
016b84c482 Mart*0421        OPPSconvectCount(I,J,K) = OPPSconvectCount(I,J,K) + 1. _d 0
6580feb7eb Jean*0422
893df04db0 Mart*0423 C     jump here if k = maxdepth or if level not unstable, go to next
                0424 C     profile point
6580feb7eb Jean*0425
893df04db0 Mart*0426  1000  continue
6580feb7eb Jean*0427
893df04db0 Mart*0428 C     end  of k-loop
6580feb7eb Jean*0429
3daafce20b Jean*0430       ENDDO
893df04db0 Mart*0431
94a46dfe0d Jean*0432 C--   End IF (DIFFERENT_MULTIPLE)
893df04db0 Mart*0433       ENDIF
                0434
                0435       RETURN
                0436       END
6580feb7eb Jean*0437
                0438 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0439
                0440       _RL FUNCTION STATE1( sLoc, tLoc, i, j, kRef, bi, bj, myThid )
                0441
893df04db0 Mart*0442 C     !DESCRIPTION: \bv
                0443 C     *===============================================================*
                0444 C     | o SUBROUTINE STATE1
                0445 C     |   Calculates rho(S,T,p)
                0446 C     |   It is absolutely necessary to compute
                0447 C     |   the full rho and not sigma=rho-rhoConst, because
                0448 C     |   density is used as a scale factor for fluxes and velocities
                0449 C     *===============================================================*
                0450 C     \ev
                0451
                0452 C     !USES:
                0453       IMPLICIT NONE
                0454 C     == Global variables ==
                0455 #include "SIZE.h"
                0456 #include "EEPARAMS.h"
                0457 #include "PARAMS.h"
                0458 #include "GRID.h"
                0459 #include "DYNVARS.h"
1a15b18a09 Jean*0460 #ifdef ALLOW_NONHYDROSTATIC
                0461 # include "NH_VARS.h"
                0462 #endif /* ALLOW_NONHYDROSTATIC */
893df04db0 Mart*0463
                0464 C     !INPUT/OUTPUT PARAMETERS:
                0465 C     == Routine arguments ==
6580feb7eb Jean*0466       INTEGER i, j, kRef, bi, bj, myThid
                0467       _RL tLoc, sLoc
893df04db0 Mart*0468
                0469 C     !LOCAL VARIABLES:
                0470 C     == Local variables ==
6580feb7eb Jean*0471       _RL rhoLoc
893df04db0 Mart*0472       _RL pLoc
                0473
                0474 CMLC     estimate pressure from depth at cell centers
                0475 CML      mtoSI = gravity*rhoConst
                0476 CML      pLoc = ABS(rC(kRef))*mtoSI
                0477
6580feb7eb Jean*0478       IF ( usingZCoords ) THEN
893df04db0 Mart*0479 C     in Z coordinates the pressure is rho0 * (hydrostatic) Potential
1a15b18a09 Jean*0480 #ifdef ALLOW_NONHYDROSTATIC
                0481        IF ( selectP_inEOS_Zc.EQ.3 ) THEN
                0482         pLoc = rhoConst*( totPhiHyd(i,j,kRef,bi,bj)
                0483      &                  + phi_nh(i,j,kRef,bi,bj)
                0484      &                  + phiRef(2*kRef)
                0485      &                  )*maskC(i,j,kRef,bi,bj)
                0486        ELSEIF ( selectP_inEOS_Zc.EQ.2 ) THEN
                0487 #else /* ALLOW_NONHYDROSTATIC */
                0488        IF     ( selectP_inEOS_Zc.EQ.2 ) THEN
                0489 #endif /* ALLOW_NONHYDROSTATIC */
893df04db0 Mart*0490 C----------
                0491 C     NOTE: For now, totPhiHyd only contains the Potential anomaly
1a15b18a09 Jean*0492 C           since PhiRef has not (yet) been added in S/R DIAGS_PHI_HYD
893df04db0 Mart*0493 C----------
                0494         pLoc = rhoConst*( totPhiHyd(i,j,kRef,bi,bj)
1a15b18a09 Jean*0495      &                  + phiRef(2*kRef)
                0496      &                  )*maskC(i,j,kRef,bi,bj)
                0497 c      ELSEIF ( selectP_inEOS_Zc.EQ.1 ) THEN
                0498 C note: for the case selectP_inEOS_Zc=0, also use pRef4EOS (set to
                0499 C       rhoConst*phiRef(2*k) ) to reproduce same previous machine truncation
                0500        ELSEIF ( selectP_inEOS_Zc.LE.1 ) THEN
                0501         pLoc = pRef4EOS(kRef)*maskC(i,j,kRef,bi,bj)
893df04db0 Mart*0502        ELSE
1a15b18a09 Jean*0503         pLoc = rhoConst*phiRef(2*kRef)*maskC(i,j,kRef,bi,bj)
893df04db0 Mart*0504        ENDIF
6580feb7eb Jean*0505       ELSEIF ( usingPCoords ) THEN
1a15b18a09 Jean*0506 C     in P coordinates the pressure is just the coordinate of the tracer point
893df04db0 Mart*0507        pLoc = rC(kRef)* maskC(i,j,kRef,bi,bj)
                0508       ENDIF
                0509
6580feb7eb Jean*0510       CALL FIND_RHO_SCALAR( tLoc, sLoc, pLoc, rhoLoc, myThid )
                0511       STATE1 = rhoLoc
893df04db0 Mart*0512
                0513 #endif /* ALLOW_OPPS */
                0514       RETURN
                0515       END
                0516
                0517 #ifdef OPPS_ORGCODE
6580feb7eb Jean*0518 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0519
                0520 C Listed below is the subroutine for use in parallel 3-d circulation code.
                0521 C It has been used in the parallel semtner-chervin code and is now being used
                0522 C In the POP code.  The subroutine is called nlopps (long story to explain why).
893df04db0 Mart*0523
6580feb7eb Jean*0524 C I've attached the version of lopps that we've been using in the simulations.
                0525 C There is one common block that is different from the standard model commons
                0526 C (countc) and it is not needed if the array convadj is not used.  The routine
                0527 C does need "kmp" which is why the boundc common is included. For massively
                0528 C parallel codes (like POP) we think this will work well when converted from a
29188f9691 Jean*0529 C "slab" (i=is,ie) to a column, which just means removing the "do i=is,ie" loop.
                0530 C There are differences between this
                0531 C code and the 1-d code and the earlier scheme implemented in 3-d models. These
                0532 C differences are described below.
893df04db0 Mart*0533
                0534       subroutine nlopps(j,is,ie,ta,sa,gcmdz)
6580feb7eb Jean*0535
893df04db0 Mart*0536       parameter (imt = 361 , jmt = 301 , km = 30 )
6580feb7eb Jean*0537
                0538 C     Nlopps:   E. Skyllingstad and T. Paluszkiewicz
                0539
                0540 C     Version: December 11, 1996
                0541
                0542 C     Nlopps:  This version of lopps is significantly different from
                0543 C     the original code developed by R. Romea and T. Paluskiewicz.  The
                0544 C     code uses a flux constraint to control the change in T and S at
                0545 C     each grid level.  First, a plume profile of T,S, and W are
                0546 C     determined using the standard plume model, but with a detraining
                0547 C     mass instead of entraining.  Thus, the T and S plume
                0548 C     characteristics still change, but the plume contracts in size
                0549 C     rather than expanding ala classical entraining plumes.  This
                0550 C     is heuristically more in line with large eddy simulation results.
                0551 C     At each grid level, the convergence of plume velocity determines
                0552 C     the flux of T and S, which is conserved by using an upstream
                0553 C     advection.  The vertical velocity is balanced so that the area
                0554 C     weighted upward velocity equals the area weighted downdraft
                0555 C     velocity, ensuring mass conservation. The present implementation
                0556 C     adjusts the plume for a time period equal to the time for 1/2 of
                0557 C     the mass of the fastest moving level to move downward.  As a
                0558 C     consequence, the model does not completely adjust the profile at
                0559 C     each model time step, but provides a smooth adjustment over time.
893df04db0 Mart*0560
                0561 c      include "params.h"
                0562 c      include "plume_fast_inc.h"
                0563 c      include "plume_fast.h"
                0564 c #include "loppsd.h"
                0565
                0566       real ta(imt,km),sa(imt,km),gcmdz(km),dz(km)
                0567       real pdensity,wsqr,oldflux,newflux,entrainrate,adtemp
                0568       REAL Del,D,dza1,dza2,kd,kd1,Smix,Thmix,PlumeS,PlumeT,PlumeD
                0569
                0570       INTEGER i,j,k
6580feb7eb Jean*0571 Clfh
893df04db0 Mart*0572       integer is,ie,k2
6580feb7eb Jean*0573 Clfh
893df04db0 Mart*0574       REAL D1,D2,state1,Density
                0575       REAL dz1,dz2
                0576       REAL radius,StartingFlux
                0577       real ttemp(km),stemp(km),taa(km),saa(km)
                0578       real wda(km),tda(km),sda(km),mda(km)
                0579       real dtts,dt,sumo,sumn
                0580       integer ntime,nn,kmx,ic
6580feb7eb Jean*0581
893df04db0 Mart*0582       LOGICAL debug,done
                0583       INTEGER MAX_ABE_ITERATIONS
                0584       PARAMETER(MAX_ABE_ITERATIONS=1)
                0585       REAL PlumeRadius
                0586       REAL STABILITY_THRESHOLD
                0587       REAL FRACTIONAL_AREA
                0588       REAL MAX_FRACTIONAL_AREA
                0589       REAL VERTICAL_VELOCITY
                0590       REAL ENTRAINMENT_RATE
                0591       REAL e2
                0592       PARAMETER ( PlumeRadius          =  100.D0   )
                0593       PARAMETER ( STABILITY_THRESHOLD  =  -1.E-4   )
                0594       PARAMETER ( FRACTIONAL_AREA      =  .1E0    )
                0595       PARAMETER ( MAX_FRACTIONAL_AREA  =  .8E0     )
                0596       PARAMETER ( VERTICAL_VELOCITY    =  .02E0   )
                0597       PARAMETER ( ENTRAINMENT_RATE     =  -.05E0     )
                0598       PARAMETER ( e2    =   2.E0*ENTRAINMENT_RATE  )
                0599       ! Arrays.
                0600       REAL Ad(km),Sd(km),Td(km),Wd(km),Dd(km),Md(km)
                0601       REAL Se(km),Te(km),We(km),De(km)
                0602       REAL PlumeEntrainment(km)
                0603       REAL GridThickness(km)
                0604
6580feb7eb Jean*0605 C input kmp through a common block
                0606
893df04db0 Mart*0607       common / boundc / wsx(imt,jmt),wsy(imt,jmt),hfs(imt,jmt),
                0608      1                  ple(imt,jmt),kmp(imt,jmt),kmq(imt,jmt)
                0609 cwmseas
                0610      &                 ,wsx1(imt,jmt),wsy1(imt,jmt)
                0611      1                 ,wsx2(imt,jmt),wsy2(imt,jmt)
6580feb7eb Jean*0612
                0613 C input the variables through a common
                0614
893df04db0 Mart*0615       logical problem
                0616       common /countc/ convadj(imt,jmt,km),ics,depth(km),problem
                0617
6580feb7eb Jean*0618 C-----may want to setup an option to get this only on first call
                0619 C     otherwise it is repetive
                0620 C     griddz is initialize by call to setupgrid
893df04db0 Mart*0621
                0622         dtts = 2400
6580feb7eb Jean*0623
893df04db0 Mart*0624         do k=1,km
                0625           dz(k) = 0.01*gcmdz(k)
                0626         enddo
6580feb7eb Jean*0627
893df04db0 Mart*0628         do k=1,km
                0629            GridThickness(k) = dz(k)
                0630         enddo
6580feb7eb Jean*0631
                0632 C modified to loop over slab
                0633
893df04db0 Mart*0634       DO i=is,ie
                0635
                0636         numgridpoints=kmp(i,j)
6580feb7eb Jean*0637
                0638 C  go to next column if only 1 grid point or on land
                0639
893df04db0 Mart*0640         if(numgridpoints.le.1) goto 1100
6580feb7eb Jean*0641
                0642 C loop over depth
                0643
893df04db0 Mart*0644       debug = .false.
6580feb7eb Jean*0645
                0646 C first save copy of initial profile
                0647
893df04db0 Mart*0648       DO k=1,NumGridPoints
                0649          stemp(k)=sa(i,k)
                0650          ttemp(k)=ta(i,k)
6580feb7eb Jean*0651
                0652 C do a check of t and s range, if out of bounds set flag
                0653
893df04db0 Mart*0654          if(problem) then
                0655             write(0,*)"Code in trouble before this nlopps call"
                0656             return
                0657          endif
6580feb7eb Jean*0658
893df04db0 Mart*0659          if(sa(i,k).gt.40..or.ta(i,k).lt.-4.0) then
                0660             problem = .true.
                0661             write(0,*)"t out of range at j ",j
                0662             debug = .true.
                0663             return
                0664          endif
                0665       ENDDO
                0666
                0667       if(debug) then
                0668         write(*,*)"T and S Profile at  ",i,j
                0669         write(*,*)(ta(i,k),sa(i,k),k=1,NumGridPoints)
                0670       endif
                0671
                0672       DO k=1,NumGridPoints-1
6580feb7eb Jean*0673
                0674 C initialize the plume T,S,density, and w velocity
                0675
893df04db0 Mart*0676           Sd(k)=stemp(k)
                0677           Td(k)=ttemp(k)
                0678           Dd(k)=state1(stemp(k),ttemp(k),k)
                0679           De(k)=Dd(k)
                0680 c          Wd(k)=VERTICAL_VELOCITY
6580feb7eb Jean*0681
                0682 C guess at initial top grid cell vertical velocity
                0683
893df04db0 Mart*0684           Wd(k) = 0.03
6580feb7eb Jean*0685
                0686 C these estimates of initial plume velocity based on plume size and
                0687 C top grid cell water mass
                0688
893df04db0 Mart*0689 c          Wd(k) = 0.5*dz(k)/(dtts*FRACTIONAL_AREA)
                0690 c          Wd(k) = 0.5*dz(k)/dtts
6580feb7eb Jean*0691
893df04db0 Mart*0692           wsqr=Wd(k)*Wd(k)
                0693           PlumeEntrainment(k) = 0.0
6580feb7eb Jean*0694
893df04db0 Mart*0695           if(debug) write(0,*) 'Doing old lowerparcel'
                0696           radius=PlumeRadius
                0697           StartingFlux=radius*radius*Wd(k)*Dd(k)
                0698           oldflux=StartingFlux
3daafce20b Jean*0699
893df04db0 Mart*0700           dz2=GridThickness(k)
                0701           DO k2=k,NumGridPoints-1
                0702             D1=state1(Sd(k2),Td(k2),k2+1)
                0703             D2=state1(stemp(k2+1),ttemp(k2+1),k2+1)
                0704             De(k2+1)=D2
6580feb7eb Jean*0705
                0706 C To start downward, parcel has to initially be heavier than environment
                0707 C but after it has started moving, we continue plume until plume tke or
                0708 C flux goes negative
                0709
893df04db0 Mart*0710             IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN
                0711                  dz1=dz2
                0712                  dz2=GridThickness(k2+1)
6580feb7eb Jean*0713
                0714 C define mass flux according to eq. 4 from paper
                0715
893df04db0 Mart*0716                  newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*0.50*
                0717      .              (dz1+dz2)
6580feb7eb Jean*0718
893df04db0 Mart*0719                  PlumeEntrainment(k2+1) = newflux/StartingFlux
6580feb7eb Jean*0720
893df04db0 Mart*0721                  IF(newflux.LT.0.0) then
                0722                      if(debug) then
                0723                       write(0,*)"Plume entrained to zero at ",k2
                0724                      endif
                0725                      maxdepth = k2
                0726                      if(maxdepth.eq.k) goto 1000
                0727                      goto 1
                0728                  endif
6580feb7eb Jean*0729
                0730 C entrainment rate is basically a scaled mass flux dM/M
                0731
893df04db0 Mart*0732                  entrainrate = (newflux - oldflux)/newflux
                0733                  oldflux = newflux
6580feb7eb Jean*0734
                0735 C mix var(s) are the average environmental values over the two grid levels
                0736
893df04db0 Mart*0737                  smix=(dz1*stemp(k2)+dz2*stemp(k2+1))/(dz1+dz2)
                0738                  thmix=(dz1*ttemp(k2)+dz2*ttemp(k2+1))/(dz1+dz2)
6580feb7eb Jean*0739
                0740 C first compute the new salinity and temperature for this level
                0741 C using equations 3.6 and 3.7 from the paper
                0742
893df04db0 Mart*0743                   sd(k2+1)=sd(k2) - entrainrate*(smix - sd(k2))
                0744                   td(k2+1)=td(k2) - entrainrate*(thmix - td(k2))
6580feb7eb Jean*0745
                0746 C compute the density at this level for the buoyancy term in the
                0747 C vertical k.e. equation
                0748
893df04db0 Mart*0749                  Dd(k2+1)=state1(Sd(k2+1),Td(k2+1),k2+1)
6580feb7eb Jean*0750
                0751 C next, solve for the vertical velocity k.e. using combined eq. 4
                0752 C and eq 5 from the paper
                0753
893df04db0 Mart*0754                  if(debug) then
a6de755a30 Jean*0755                   write(0,*)"Dd,De,entr,k ",Dd(k2),De(k2),entrainrate,k2
893df04db0 Mart*0756                  endif
6580feb7eb Jean*0757
893df04db0 Mart*0758                  wsqr = wsqr - wsqr*abs(entrainrate)+ 9.81*
                0759      .             (dz1*(Dd(k2)-De(k2))/De(k2)
                0760      .             +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1))
6580feb7eb Jean*0761
                0762 C if negative k.e. then plume has reached max depth, get out of loop
                0763
893df04db0 Mart*0764                  IF(wsqr.LT.0.0)then
                0765                      maxdepth = k2
                0766                      if(debug) then
                0767                       write(0,*)"Plume velocity went to zero at ",k2
                0768                      endif
                0769                      if(maxdepth.eq.k) goto 1000
                0770                      goto 1
                0771                  endif
                0772                  Wd(k2+1)=sqrt(wsqr)
6580feb7eb Jean*0773
                0774 C compute a new radius based on the new mass flux at this grid level
                0775
893df04db0 Mart*0776                  radius=sqrt(newflux/(Wd(k2)*Dd(k2)))
                0777               ELSE
                0778                  maxdepth=k2
                0779                  if(maxdepth.eq.k) goto 1000
                0780                  GOTO 1
                0781               ENDIF
                0782           ENDDO
6580feb7eb Jean*0783
                0784 C plume has reached the bottom
                0785
893df04db0 Mart*0786           MaxDepth=NumGridPoints
6580feb7eb Jean*0787
893df04db0 Mart*0788 1         continue
6580feb7eb Jean*0789
893df04db0 Mart*0790           Ad(k)=FRACTIONAL_AREA
                0791           IC=0
6580feb7eb Jean*0792
                0793 C start iteration on fractional area, not used in OGCM implementation
                0794
893df04db0 Mart*0795           DO IC=1,Max_ABE_Iterations
6580feb7eb Jean*0796
                0797 C next compute the mass flux beteen each grid box using the entrainment
                0798
893df04db0 Mart*0799  92          continue
                0800              Md(k)=Wd(k)*Ad(k)
6580feb7eb Jean*0801
893df04db0 Mart*0802              DO k2=k+1,MaxDepth
                0803                Md(k2)=Md(k)*PlumeEntrainment(k2)
                0804                if(debug) then
                0805                  write(0,*)"Md, Wd, and  k are ",Md(k2),Wd(k2),k2
                0806                endif
                0807              ENDDO
6580feb7eb Jean*0808
                0809 C Now move on to calculate new temperature using flux from
                0810 C Td, Sd, Wd, ta, sa, and we. Values for these variables are at
                0811 C center of grid cell, use weighted average to get boundary values
                0812
                0813 C use a timestep limited by the GCM model timestep and the maximum plume
                0814 C velocity (CFL criteria)
                0815
                0816 C calculate the weighted wd, td, and sd
                0817
893df04db0 Mart*0818              dt = dtts
                0819              do k2=k,maxdepth-1
                0820                 dt = min(dt,dz(k2)/wd(k2))
6580feb7eb Jean*0821
                0822 C time integration will be integer number of steps to get one
                0823 C gcm time step
                0824
893df04db0 Mart*0825                 ntime = nint(0.5*int(dtts/dt))
                0826                 if(ntime.eq.0) then
                0827                    ntime = 1
                0828                 endif
6580feb7eb Jean*0829
                0830 C make sure area weighted vertical velocities match; in other words
                0831 C make sure mass in equals mass out at the intersection of each grid
                0832 C cell.
                0833
893df04db0 Mart*0834                 mda(k2) = (md(k2)*dz(k2)+md(k2+1)*dz(k2+1))/
                0835      *                    (dz(k2)+dz(k2+1))
6580feb7eb Jean*0836
893df04db0 Mart*0837                 wda(k2) = (wd(k2)*dz(k2)+wd(k2+1)*dz(k2+1))/
                0838      *                    (dz(k2)+dz(k2+1))
6580feb7eb Jean*0839
893df04db0 Mart*0840                 tda(k2) = td(k2)
                0841                 sda(k2) = sd(k2)
6580feb7eb Jean*0842
893df04db0 Mart*0843                 taa(k2) = ttemp(k2+1)
                0844                 saa(k2) = stemp(k2+1)
6580feb7eb Jean*0845
893df04db0 Mart*0846              enddo
                0847              dt = min(dt,dtts)
                0848              if(debug) then
                0849                write(0,*)"Time step is ", dt
                0850              endif
                0851              tda(maxdepth) = td(maxdepth)
                0852              sda(maxdepth) = sd(maxdepth)
6580feb7eb Jean*0853
                0854 C do top and bottom points first
                0855
893df04db0 Mart*0856              kmx = maxdepth-1
                0857              do nn=1,ntime
                0858
                0859                ttemp(k) =  ttemp(k)-
                0860      *                  (mda(k)*(tda(k)-taa(k)))*dt*recip_drF(k)
6580feb7eb Jean*0861
893df04db0 Mart*0862                stemp(k) =  stemp(k)-
                0863      *                  (mda(k)*(sda(k)-saa(k)))*dt*recip_drF(k)
6580feb7eb Jean*0864
                0865 C now do inner points if there are any
                0866
893df04db0 Mart*0867                if(Maxdepth-k.gt.1) then
                0868                  do k2=k+1,Maxdepth-1
6580feb7eb Jean*0869
893df04db0 Mart*0870                    ttemp(k2) = ttemp(k2) +
                0871      *              (mda(k2-1)*(tda(k2-1)-taa(k2-1))-
                0872      *              mda(k2)*(tda(k2)-taa(k2)))
                0873      *              *dt*recip_drF(k2)
                0874
                0875                   stemp(k2) = stemp(k2) +
                0876      *              (mda(k2-1)*(sda(k2-1)-saa(k2-1))-
                0877      *              mda(k2)*(sda(k2)-saa(k2)))
                0878      *              *dt*recip_drF(k2)
                0879
                0880                  enddo
                0881                endif
                0882                ttemp(kmx+1) =  ttemp(kmx+1)+
                0883      *                  (mda(kmx)*(tda(kmx)-taa(kmx)))*
                0884      *                  dt*recip_drF(kmx+1)
6580feb7eb Jean*0885
893df04db0 Mart*0886                stemp(kmx+1) =  stemp(kmx+1)+
                0887      *                  (mda(kmx)*(sda(kmx)-saa(kmx)))*
                0888      *                  dt*recip_drF(kmx+1)
6580feb7eb Jean*0889
                0890 C set the environmental temp and salinity to equal new fields
                0891
893df04db0 Mart*0892                 do k2=1,maxdepth-1
                0893                   taa(k2) = ttemp(k2+1)
                0894                   saa(k2) = stemp(k2+1)
                0895                 enddo
6580feb7eb Jean*0896
                0897 C end loop on number of time integration steps
                0898
893df04db0 Mart*0899              enddo
                0900           ENDDO
                0901 999       continue
6580feb7eb Jean*0902
                0903 C assume that it converged, so update the ta and sa with new fields
                0904
893df04db0 Mart*0905 c          if(i.gt.180.and.j.gt.200.and.i.lt.240) then
                0906 c            write(*,*)"Converged ",i,j,k,maxdepth,ttemp(k+1),ta(i,k+1)
                0907 c          endif
                0908           do k2=k,maxdepth
                0909             convadj(i,j,k2) = convadj(i,j,k2) + (ttemp(k2)-
                0910      *          ta(i,k2))
                0911             sa(i,k2) = stemp(k2)
                0912             ta(i,k2) = ttemp(k2)
                0913 c          if(i.gt.180.and.j.gt.200.and.i.lt.240) then
                0914 c            write(*,*)"convadj ",convadj(i,j,k2)
                0915 c          endif
6580feb7eb Jean*0916
                0917 C see if nlopps messed up
                0918
893df04db0 Mart*0919             if(sa(i,k).gt.40..or.ta(i,k).lt.-4.0) then
                0920                problem = .true.
                0921                write(0,*)"t out of range at j after adjust",j
                0922                debug = .true.
                0923             endif
6580feb7eb Jean*0924
893df04db0 Mart*0925           enddo
6580feb7eb Jean*0926
                0927 C jump here if k = maxdepth or if level not unstable, go to next
                0928 C profile point
                0929
893df04db0 Mart*0930 1000      continue
6580feb7eb Jean*0931
                0932 C end loop on k, move on to next possible plume
                0933
3daafce20b Jean*0934       ENDDO
893df04db0 Mart*0935 1100  continue
6580feb7eb Jean*0936
                0937 C i loop
                0938
893df04db0 Mart*0939       ENDDO
                0940       END
                0941
6611306112 Ed H*0942 #endif /* OPPS_ORGCODE */