pkg/salt_plume/salt_plume_volfrac.F

e10b4e6008 Patr*0001 #include "SALT_PLUME_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SALT_PLUME_VOLFRAC
                0005 C     !INTERFACE:
                0006       SUBROUTINE SALT_PLUME_VOLFRAC(
                0007      I                       bi, bj, myTime, myIter, myThid )
                0008
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | SUBROUTINE SALT_PLUME_VOLFRAC
                0012 C     | o Compute saltplume penetration.
                0013 C     *==========================================================*
ba0b047096 Mart*0014 C     | Compute fraction of volume flux associated with saltplume
                0015 C     | flux penetrating through the entire water columns due to
e10b4e6008 Patr*0016 C     | rejected salt during freezing.
ba0b047096 Mart*0017 C     |
e10b4e6008 Patr*0018 C     | For example, if surface value is Saltplume0,
                0019 C     | and each level gets equal fraction 1/5 down to SPDepth=5,
ba0b047096 Mart*0020 C     | SALT_PLUME_VOLFRAC will report
e10b4e6008 Patr*0021 C     | dSPvolkLev2Above[2to1,3to2,4to3,5to4,6to5] = [4/5,3/5,2/5,1/5,  0]
                0022 C     | dSPvolSurf2kLev [1to1,1to2,1to3,1to4,1to5] = [1/5,1/5,1/5,1/5,1/5]
                0023 C     | sum [into5] = 1to5 + 6to5 - 5to4 = 1/5 +   0 - 1/5 = 0
                0024 C     |     [into4] = 1to4 + 5to4 - 4to3 = 1/5 + 1/5 - 2/5 = 0
                0025 C     |     [into3] = 1to3 + 4to3 - 3to2 = 1/5 + 2/5 - 3/5 = 0
                0026 C     |     [into2] = 1to2 + 3to2 - 2to1 = 1/5 + 3/5 - 4/5 = 0
                0027 C     |     [into1] = 1to1 + 2to1 - 1to[1,2,3,4,5] = 1/5 + 4/5 - 5/5 = 0
                0028 C     | NOTE: volume will always be conserved.
                0029 C     | =====
                0030 C     | Written by   : ATN (based on SALT_PLUME_FRAC)
                0031 C     | Date         : Apr 14, 2014
                0032 C     *==========================================================*
                0033 C     \ev
                0034
                0035 C     !USES:
                0036       IMPLICIT NONE
                0037 #include "SIZE.h"
                0038 #include "GRID.h"
                0039 #include "SALT_PLUME.h"
                0040 #include "EEPARAMS.h"
                0041 #include "PARAMS.h"
                0042
                0043 C     !INPUT/OUTPUT PARAMETERS:
                0044 C     input arguments
                0045 C     SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point
                0046 C     myTime  :: Current time in simulation
                0047 C     myIter  :: Current iteration number in simulation
                0048 C     myThid  :: My Thread Id. number
                0049       INTEGER bi,bj
                0050       _RL     myTime
                0051       INTEGER myIter
                0052       INTEGER myThid
                0053 C     input/output arguments
                0054 C      CHARACTER*(MAX_LEN_MBUF) msgBuf
                0055 CEOP
                0056
                0057 #ifdef ALLOW_SALT_PLUME
                0058 #ifdef SALT_PLUME_VOLUME
                0059
                0060 C     !LOCAL VARIABLES:
                0061       _RL     dMbdt        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0062       _RL     temp         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0063       _RL     dplumek
                0064       INTEGER SPkBottom    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0065       INTEGER i,j,k,kp1,Nlev,Nrp1
                0066       INTEGER imt
                0067       parameter( imt=(sNx+2*OLx)*(sNy+2*OLy) )
                0068
                0069 C initialize at every time step
                0070       Nrp1=Nr+1
                0071       DO k=1,Nr
                0072        DO j=1-OLy,sNy+OLy
                0073         DO i=1-OLx,sNx+OLx
                0074          dSPvolSurf2kLev(i,j,k,bi,bj)  = 0. _d 0
                0075          dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
                0076          SPplumek(i,j,k,bi,bj)         = 1. _d 0
                0077         ENDDO
                0078        ENDDO
                0079       ENDDO
                0080       DO j=1-OLy,sNy+OLy
                0081        DO i=1-OLx,sNx+OLx
                0082         SPplumek(i,j,Nrp1,bi,bj)         = 1. _d 0
                0083         SPbrineVolFlux(i,j,bi,bj)        = 0. _d 0
                0084         SPkBottom(i,j)                   = 0
                0085        ENDDO
                0086       ENDDO
                0087
                0088 C call salt_plume_frac to fill in SPplumek and SPkBottom
                0089 C use dMbdt+temp as a temporary arrays here to save memory:
                0090       DO k = Nrp1,1,-1
                0091        DO j=1-Oly,sNy+Oly
                0092         DO i=1-Olx,sNx+Olx
                0093          temp(i,j)=SaltPlumeDepth(i,j,bi,bj)
                0094          dMbdt(i,j)=abs(rF(k))
                0095         ENDDO
                0096        ENDDO
                0097        CALL SALT_PLUME_FRAC(
                0098      I               imt,oneRS,temp,
                0099 #ifdef SALT_PLUME_SPLIT_BASIN
d74df16944 An T*0100      I               XC(1-Olx,1-Oly,bi,bj),YC(1-Olx,1-Oly,bi,bj),
e10b4e6008 Patr*0101 #endif
                0102      U               dMbdt,
                0103      I               myTime, 1, myThid )
                0104        DO j=1-Oly,sNy+Oly
                0105         DO i=1-Olx,sNx+Olx
                0106          SPplumek(i,j,k,bi,bj)=dMbdt(i,j)
                0107          IF(SPplumek(i,j,k,bi,bj).GT. 0.9999999) THEN
                0108           SPkBottom(i,j)=k
                0109          ENDIF
                0110         ENDDO
                0111        ENDDO
                0112       ENDDO
                0113
                0114 C reinitialize dMbdt = 0
                0115       DO j=1-Oly,sNy+Oly
                0116        DO i=1-Olx,sNx+Olx
                0117         dMbdt(i,j)=0. _d 0
                0118        ENDDO
                0119       ENDDO
                0120
                0121 C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev
                0122 C units:
ba0b047096 Mart*0123 C Sbrine=dsb/dt*dt/(rhoConst*SPalpha*drF)[g/m2/s*s/(kg/m3*m)]=[g/kg]
e10b4e6008 Patr*0124 C SPplumek : fraction : unitless
ba0b047096 Mart*0125 C SaltPlumeFlux: dsb/dt [g/kg.kg/m^2/s = g/m^2/s]
e10b4e6008 Patr*0126 C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]
                0127 C                        = dsb/dt / (dsb/dt*dt/(rhoConst*SPalpha*drF))
                0128 C                        = rhoConst*SPalpha*drF/dt [kg/m3 m/s]=[kg/m2/s]
                0129 C dVbrine/dt = dMb/dt 1/rhoConst [m/s]
                0130
ba0b047096 Mart*0131 C has 2 ways to define brine properties: either provide
e10b4e6008 Patr*0132 C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.
ba0b047096 Mart*0133 C (A) SPalpha:  can calc SPbrineSalt as fxn of dhice/dt,
e10b4e6008 Patr*0134 C     constrained by SPbrineSaltmax:
                0135 C     SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF(1)*dt
                0136 C     SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)
                0137 C     dMbdt = saltPlumeFlux / SPbrineSalt
                0138 C           = rhoConst*SPalpha*drF(1)/dt <-- a function of SPalpha
ba0b047096 Mart*0139 C (B) SPbrinesalt provided
e10b4e6008 Patr*0140 C     dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt
                0141
                0142 C Assuming we go with (B) here:
                0143       DO j=1-OLy,sNy+OLy
                0144        DO i=1-OLx,sNx+OLx
                0145 C brine mass and volume at surface:
                0146         dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSconst
                0147         SPbrineVolFlux(i,j,bi,bj)=dMbdt(i,j)*mass2rUnit
                0148        ENDDO
                0149       ENDDO
                0150
                0151 C Distributing down: this is always from level 1 to depth
                0152       DO k=Nr,1,-1
                0153        DO j=1-OLy,sNy+OLy
                0154         DO i=1-OLx,sNx+OLx
                0155          dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
                0156          dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*SPbrineVolFlux(i,j,bi,bj)
                0157         ENDDO
                0158        ENDDO
                0159       ENDDO
                0160
                0161 C Now volume up: need to scan from bottom of SPDepth
                0162       DO j=1-OLy,sNy+OLy
                0163        DO i=1-OLx,sNx+OLx
                0164         Nlev=SPkBottom(i,j)
                0165         IF(Nlev.GE.1 .AND. Nlev.LE.Nr) THEN
                0166          DO k=Nlev,1,-1
                0167           kp1=k+1
                0168           dSPvolkLev2Above(i,j,k,bi,bj)=dSPvolkLev2Above(i,j,kp1,bi,bj)
                0169      &                                  - dSPvolSurf2kLev(i,j,k,bi,bj)
                0170          ENDDO
                0171         ENDIF
                0172        ENDDO
                0173       ENDDO
                0174
                0175 #ifdef ALLOW_DIAGNOSTICS
                0176       IF ( useDiagnostics ) THEN
                0177        CALL DIAGNOSTICS_FILL(
                0178      &      SPplumek,'PLUMEKB1',0,Nr,1,bi,bj,myThid )
                0179       ENDIF
                0180 #endif /* ALLOW_DIAGNOSTICS */
                0181
                0182 #endif /* SALT_PLUME_VOLUME */
                0183 #endif /* ALLOW_SALT_PLUME */
                0184
                0185       RETURN
                0186       END