C $Header: /u/gcmpack/MITgcm/pkg/salt_plume/salt_plume_volfrac.F,v 1.2 2014/05/21 13:59:39 atn Exp $
C $Name: $
#include "SALT_PLUME_OPTIONS.h"
CBOP
C !ROUTINE: SALT_PLUME_VOLFRAC
C !INTERFACE:
SUBROUTINE SALT_PLUME_VOLFRAC(
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE SALT_PLUME_VOLFRAC
C | o Compute saltplume penetration.
C *==========================================================*
C | Compute fraction of volume flux associated with saltplume
C | flux penetrating through the entire water columns due to
C | rejected salt during freezing.
C |
C | For example, if surface value is Saltplume0,
C | and each level gets equal fraction 1/5 down to SPDepth=5,
C | SALT_PLUME_VOLFRAC will report
C | dSPvolkLev2Above[2to1,3to2,4to3,5to4,6to5] = [4/5,3/5,2/5,1/5, 0]
C | dSPvolSurf2kLev [1to1,1to2,1to3,1to4,1to5] = [1/5,1/5,1/5,1/5,1/5]
C | sum [into5] = 1to5 + 6to5 - 5to4 = 1/5 + 0 - 1/5 = 0
C | [into4] = 1to4 + 5to4 - 4to3 = 1/5 + 1/5 - 2/5 = 0
C | [into3] = 1to3 + 4to3 - 3to2 = 1/5 + 2/5 - 3/5 = 0
C | [into2] = 1to2 + 3to2 - 2to1 = 1/5 + 3/5 - 4/5 = 0
C | [into1] = 1to1 + 2to1 - 1to[1,2,3,4,5] = 1/5 + 4/5 - 5/5 = 0
C | NOTE: volume will always be conserved.
C | =====
C | Written by : ATN (based on SALT_PLUME_FRAC)
C | Date : Apr 14, 2014
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
#include "SIZE.h"
#include "GRID.h"
#include "SALT_PLUME.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
C !INPUT/OUTPUT PARAMETERS:
C input arguments
C SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: My Thread Id. number
INTEGER bi,bj
_RL myTime
INTEGER myIter
INTEGER myThid
C input/output arguments
C CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_VOLUME
C !LOCAL VARIABLES:
_RL dMbdt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL temp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dplumek
INTEGER SPkBottom (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER i,j,k,kp1,Nlev,Nrp1
INTEGER imt
parameter( imt=(sNx+2*OLx)*(sNy+2*OLy) )
C initialize at every time step
Nrp1=Nr+1
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dSPvolSurf2kLev(i,j,k,bi,bj) = 0. _d 0
dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
SPplumek(i,j,k,bi,bj) = 1. _d 0
ENDDO
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
SPplumek(i,j,Nrp1,bi,bj) = 1. _d 0
SPbrineVolFlux(i,j,bi,bj) = 0. _d 0
SPkBottom(i,j) = 0
ENDDO
ENDDO
C call salt_plume_frac to fill in SPplumek and SPkBottom
C use dMbdt+temp as a temporary arrays here to save memory:
DO k = Nrp1,1,-1
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
temp(i,j)=SaltPlumeDepth(i,j,bi,bj)
dMbdt(i,j)=abs(rF(k))
ENDDO
ENDDO
CALL SALT_PLUME_FRAC(
I imt,oneRS,temp,
#ifdef SALT_PLUME_SPLIT_BASIN
I XC(1-Olx,1-Oly,bi,bj),YC(1-Olx,1-Oly,bi,bj),
#endif
U dMbdt,
I myTime, 1, myThid )
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
SPplumek(i,j,k,bi,bj)=dMbdt(i,j)
IF(SPplumek(i,j,k,bi,bj).GT. 0.9999999) THEN
SPkBottom(i,j)=k
ENDIF
ENDDO
ENDDO
ENDDO
C reinitialize dMbdt = 0
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
dMbdt(i,j)=0. _d 0
ENDDO
ENDDO
C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev
C units:
C Sbrine=dsb/dt*dt/(rhoConst*SPalpha*drF)[psu kg/m2/s*s/(kg/m3*m)]=[psu]
C SPplumek : fraction : unitless
C SaltPlumeFlux: dsb/dt [psu.kg/m^2/s = g/m^2/s]
C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]
C = dsb/dt / (dsb/dt*dt/(rhoConst*SPalpha*drF))
C = rhoConst*SPalpha*drF/dt [kg/m3 m/s]=[kg/m2/s]
C dVbrine/dt = dMb/dt 1/rhoConst [m/s]
C has 2 ways to define brine properties: either provide
C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.
C (A) SPalpha: can calc SPbrineSalt as fxn of dhice/dt,
C constrained by SPbrineSaltmax:
C SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF(1)*dt
C SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)
C dMbdt = saltPlumeFlux / SPbrineSalt
C = rhoConst*SPalpha*drF(1)/dt <-- a function of SPalpha
C (B) SPbrinesalt provided
C dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt
C Assuming we go with (B) here:
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C brine mass and volume at surface:
dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSconst
SPbrineVolFlux(i,j,bi,bj)=dMbdt(i,j)*mass2rUnit
ENDDO
ENDDO
C Distributing down: this is always from level 1 to depth
DO k=Nr,1,-1
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*SPbrineVolFlux(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
C Now volume up: need to scan from bottom of SPDepth
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
Nlev=SPkBottom(i,j)
IF(Nlev.GE.1 .AND. Nlev.LE.Nr) THEN
DO k=Nlev,1,-1
kp1=k+1
dSPvolkLev2Above(i,j,k,bi,bj)=dSPvolkLev2Above(i,j,kp1,bi,bj)
& - dSPvolSurf2kLev(i,j,k,bi,bj)
ENDDO
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(
& SPplumek,'PLUMEKB1',0,Nr,1,bi,bj,myThid )
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* SALT_PLUME_VOLUME */
#endif /* ALLOW_SALT_PLUME */
RETURN
END