/[MITgcm]/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F
ViewVC logotype

Diff of /MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F

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

revision 1.1 by atn, Sun Apr 20 04:03:07 2014 UTC revision 1.4 by atn, Tue Apr 29 06:49:40 2014 UTC
# Line 7  CBOP Line 7  CBOP
7  C     !ROUTINE: SALT_PLUME_VOLFRAC  C     !ROUTINE: SALT_PLUME_VOLFRAC
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE SALT_PLUME_VOLFRAC(        SUBROUTINE SALT_PLUME_VOLFRAC(
10       I                  SaltPlumeDepth,SaltPlumeFlux,       I                       bi, bj, myTime, myIter, myThid )
      O                  SPbrineSalt,dSPvolSurf2kLev,  
      O                  dSPvolkLev2Above,dSPvolBelow2kLev,  
      I                  myTime, myIter, myThid )  
11    
12  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
13  C     *==========================================================*  C     *==========================================================*
# Line 43  C     !USES: Line 40  C     !USES:
40  #include "SIZE.h"  #include "SIZE.h"
41  #include "GRID.h"  #include "GRID.h"
42  #include "SALT_PLUME.h"  #include "SALT_PLUME.h"
43    #include "EEPARAMS.h"
44    #include "PARAMS.h"
45    
46  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
47  C     input arguments  C     input arguments
 C     imax    :: number of vertical grid points  
48  C     SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point  C     SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point
49  C     myTime  :: Current time in simulation  C     myTime  :: Current time in simulation
50  C     myIter  :: Current iteration number in simulation  C     myIter  :: Current iteration number in simulation
51  C     myThid  :: My Thread Id. number  C     myThid  :: My Thread Id. number
52        INTEGER imax        INTEGER bi,bj
53        _RL     myTime        _RL     myTime
54        INTEGER myIter        INTEGER myIter
55        INTEGER myThid        INTEGER myThid
56  C     input/output arguments  C     input/output arguments
57    C      CHARACTER*(MAX_LEN_MBUF) msgBuf
58  CEOP  CEOP
59    
60  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
61    #ifdef SALT_PLUME_VOLUME
62    
63  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
64        _RL     plumek       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)        _RL     dMbdt        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65          _RL     dVbdt        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66        _RL     dplumek        _RL     dplumek
       _RL     kLevCbot     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
67        _RL locz, zo, z20        _RL locz, zo, z20
68        INTEGER i        INTEGER i,j,k,kk,Nlev
69        _RL     one, two, three, S, So, zero        _RL     one, two, three, S, So, zero
70        parameter( one = 1. _d 0, two = 2. _d 0, three = 3. _d 0 )        parameter( one = 1. _d 0, two = 2. _d 0, three = 3. _d 0 )
71        parameter( zero = 0. _d 0 )        parameter( zero = 0. _d 0 )
# Line 73  C     This is an abbreviation of 1./(exp Line 73  C     This is an abbreviation of 1./(exp
73        _RL     recip_expOneM1        _RL     recip_expOneM1
74        parameter( recip_expOneM1 = 0.581976706869326343 )        parameter( recip_expOneM1 = 0.581976706869326343 )
75    
76  Catn: define dMbdt and dVbdt as 3d for cases of diagnostics  C initialize at every time step
       _RL dMbdt            (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL dVbdt            (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
   
   
 C initialize  
77        dplumek = 0. _d 0        dplumek = 0. _d 0
78        DO k=1,Nr        DO k=1,Nr
79         DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
80          DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
81           dSPvolSurf2kLev(i,j,k)= 0. _d 0           dSPvolSurf2kLev(i,j,k,bi,bj)  = 0. _d 0
82           dSPvolBelow2kLev(i,j,k)  = 0. _d 0           dSPvolBelow2kLev(i,j,k,bi,bj) = 0. _d 0
83           dSPvolkLev2Above(i,j,k) = 0. _d 0           dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
84           plumek(i,j,k)       = 1. _d 0           SPplumek(i,j,k,bi,bj)         = 1. _d 0
85           IF(k.EQ.Nr) THEN           IF(k.EQ.Nr) THEN
86            plumek(i,j,k+1)      = 1. _d 0            SPplumek(i,j,k+1,bi,bj)         = 1. _d 0
87            dSPvolBelow2kLev(i,j,k+1) = 0. _d 0            dSPvolBelow2kLev(i,j,k+1,bi,bj) = 0. _d 0
88           ENDIF           ENDIF
89          ENDDO          ENDDO
90         ENDDO         ENDDO
91        ENDDO        ENDDO
92        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
93         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
94          kLevCbot(i,j) = 0          SPkBottom(i,j,bi,bj)   = 0
95          SPbrineSalt(i,j) = 0. _d 0          SPbrineSalt(i,j,bi,bj) = 0. _d 0
96          dMbdt(i,j) = 0. _d 0          dMbdt(i,j)             = 0. _d 0
97          dVbdt(i,j) = 0. _d 0          dVbdt(i,j)             = 0. _d 0
98         ENDDO         ENDDO
99        ENDDO        ENDDO
100    
101        DO k = 1,Nr+1        DO k = 1,Nr+1
102         DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
103          DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
104           zloc = abs(rF(k))           locz = abs(rF(k))
105    
106           IF ( SaltPlumeDepth(i,j).GT.zloc .and.           IF ( SaltPlumeDepth(i,j,bi,bj).GT.locz .and.
107       &        SaltPlumeDepth(i,j).GT.zero ) THEN       &        SaltPlumeDepth(i,j,bi,bj).GT.zero ) THEN
108    
109            kLevCbot(i,j)=k            SPkBottom(i,j,bi,bj)=k
110  C     Default: uniform distribution, PlumeMethod=1, Npower=0  C     Default: uniform distribution, PlumeMethod=1, Npower=0
111            IF (PlumeMethod .EQ. 1) THEN            IF (PlumeMethod .EQ. 1) THEN
112             zo = abs(SaltPlumeDepth(i,j))             zo = abs(SaltPlumeDepth(i,j,bi,bj))
113             S  = one                  !input depth temp             S  = one                  !input depth temp
114             So = one             So = one
115             DO kk=1,Npower+1             DO kk=1,Npower+1
116              S  = locz*S              !raise to the Npower+1              S  = locz*S              !raise to the Npower+1
117              So = zo*So              So = zo*So
118             ENDDO             ENDDO
119             plumek(i,j,k) = max(zero,S/So)             SPplumek(i,j,k,bi,bj) = max(zero,S/So)
120    
121            ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution            ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution
122             z20 = abs(SaltPlumeDepth(i,j))             z20 = abs(SaltPlumeDepth(i,j,bi,bj))
123             S   = exp(locz/z20)-one             S   = exp(locz/z20)-one
124             So  = recip_expOneM1      !So = exp(one)-one             So  = recip_expOneM1      !So = exp(one)-one
125             plumek(i,j,k) = max(zero,S*So)             SPplumek(i,j,k,bi,bj) = max(zero,S*So)
126    
127  C     PlumeMethod = 3, distribute salt LINEARLY between SPDepth and  C     PlumeMethod = 3, distribute salt LINEARLY between SPDepth and
128  C     SPDepth/SPovershoot  C     SPDepth/SPovershoot
129  C     (1-SPovershoot)percent has already been taken into account in  C     (1-SPovershoot)percent has already been taken into account in
130  C     SPDepth calculation, i.e., SPDepth = SPovershoot*SPDepth.  C     SPDepth calculation, i.e., SPDepth = SPovershoot*SPDepth.
131            ELSEIF (PlumeMethod .EQ. 3) THEN !overshoot 20%            ELSEIF (PlumeMethod .EQ. 3) THEN !overshoot 20%
132             z20 = abs(SaltPlumeDepth(i,j))             z20 = abs(SaltPlumeDepth(i,j,bi,bj))
133             zo  = z20/SPovershoot             zo  = z20/SPovershoot
134             So  = z20-zo             So  = z20-zo
135             S   = locz-zo             S   = locz-zo
136             IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN             IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
137              plumek(i,j,k) = max(zero,S/So)              SPplumek(i,j,k,bi,bj) = max(zero,S/So)
138             ELSE             ELSE
139              plumek(i,j,k) = zero              SPplumek(i,j,k,bi,bj) = zero
140             ENDIF             ENDIF
141    
142  C     PlumeMethod = 5, dumping all salt at the top layer  C     PlumeMethod = 5, dumping all salt at the top layer
# Line 149  C     PlumeMethod = 5, dumping all salt Line 144  C     PlumeMethod = 5, dumping all salt
144             z20  = one             z20  = one
145             zo   = zero             zo   = zero
146             IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN             IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
147              plumek(i,j,k) = zero              SPplumek(i,j,k,bi,bj) = zero
148             ELSE             ELSE
149              plumek(i,j,k) = one              SPplumek(i,j,k,bi,bj) = one
150             ENDIF             ENDIF
151            ELSEIF (PlumeMethod .EQ. 6) THEN            ELSEIF (PlumeMethod .EQ. 6) THEN
152  C     PLumeMethod = 6, currently only works for Npower = 1 and 2.  C     PLumeMethod = 6, currently only works for Npower = 1 and 2.
153             z20 = abs(SaltPlumeDepth(i,j))             z20 = abs(SaltPlumeDepth(i,j,bi,bj))
154             S   = one                 !input depth temp             S   = one                 !input depth temp
155             So  = one             So  = one
156             DO kk=1,Npower+1             DO kk=1,Npower+1
# Line 163  C     PLumeMethod = 6, currently only wo Line 158  C     PLumeMethod = 6, currently only wo
158              So = z20*So              So = z20*So
159             ENDDO             ENDDO
160             IF(Npower.EQ.1) THEN      !Npower=1             IF(Npower.EQ.1) THEN      !Npower=1
161              plumek(i,j,k) = max(zero,two/z20*locz-S/So)              SPplumek(i,j,k,bi,bj) = max(zero,two/z20*locz-S/So)
162             ELSE                      !Npower=2             ELSE                      !Npower=2
163              plumek(i,j,k) = max(zero,              SPplumek(i,j,k,bi,bj) = max(zero,
164       &           three/z20*locz - three/(z20*z20)*locz*locz + S/So)       &           three/z20*locz - three/(z20*z20)*locz*locz + S/So)
165             ENDIF             ENDIF
166                        
# Line 179  C        0 -> parabola centered at SPDep Line 174  C        0 -> parabola centered at SPDep
174  C        + -> parabola offset, salt @ surface < @ SPDepth  C        + -> parabola offset, salt @ surface < @ SPDepth
175  C        - -> parabola offset, salt @ surface > @ SPDepth  C        - -> parabola offset, salt @ surface > @ SPDepth
176  C        S and So are dummy variables  C        S and So are dummy variables
177             z20 = abs(SaltPlumeDepth(i,j))             z20 = abs(SaltPlumeDepth(i,j,bi,bj))
178             zo  = z20*(one/two-Npower/200. _d 0)             zo  = z20*(one/two-Npower/200. _d 0)
179             So  = (z20*z20*z20/three)             So  = (z20*z20*z20/three)
180       &              - (z20*z20)*zo       &              - (z20*z20)*zo
# Line 187  C        S and So are dummy variables Line 182  C        S and So are dummy variables
182             S    = (locz*locz*locz/three)             S    = (locz*locz*locz/three)
183       &               - (locz*locz)*zo       &               - (locz*locz)*zo
184       &               + (locz)     *zo*zo       &               + (locz)     *zo*zo
185             plumek(i,j,k) = max(zero,(S/So))             SPplumek(i,j,k,bi,bj) = max(zero,(S/So))
186    
187            ELSE            ELSE
188             WRITE(*,*) 'salt_plume_frac: PLumeMethod =', PLumeMethod,             WRITE(*,*) 'salt_plume_frac: PLumeMethod =', PLumeMethod,
# Line 195  C        S and So are dummy variables Line 190  C        S and So are dummy variables
190             STOP 'ABNORMAL in S/R SALT_PLUME_FRAC'             STOP 'ABNORMAL in S/R SALT_PLUME_FRAC'
191            ENDIF            ENDIF
192           ELSE           ELSE
193            plumek(i,j,k) = one            SPplumek(i,j,k,bi,bj) = one
194           ENDIF           ENDIF
195          ENDDO          ENDDO
196         ENDDO         ENDDO
# Line 204  C        S and So are dummy variables Line 199  C        S and So are dummy variables
199  C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev  C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev
200  C units:  C units:
201  C Sbrine=dsb/dt*dt/(rhoConst*SPalpha*drF)[psu kg/m2/s*s/(kg/m3*m)]=[psu]  C Sbrine=dsb/dt*dt/(rhoConst*SPalpha*drF)[psu kg/m2/s*s/(kg/m3*m)]=[psu]
202  C plumekb : fraction : unitless  C SPplumek : fraction : unitless
203  C SaltPlumeFlux: dsb/dt [psu.kg/m^2/s = g/m^2/s]  C SaltPlumeFlux: dsb/dt [psu.kg/m^2/s = g/m^2/s]
204  C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]  C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]
205  C                        = dsb/dt / (dsb/dt*dt/(rhoConst*SPalpha*drF))  C                        = dsb/dt / (dsb/dt*dt/(rhoConst*SPalpha*drF))
# Line 215  C has 2 ways to define brine properties: Line 210  C has 2 ways to define brine properties:
210  C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.  C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.
211  C (A) SPalpha:  can calc SPbrineSalt as fxn of dhice/dt,  C (A) SPalpha:  can calc SPbrineSalt as fxn of dhice/dt,
212  C     constrained by SPbrineSaltmax:  C     constrained by SPbrineSaltmax:
213  C     SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF*dt  C     SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF(1)*dt
214  C     SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)  C     SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)
215  C     dMbdt = saltPlumeFlux / SPbrineSalt  C     dMbdt = saltPlumeFlux / SPbrineSalt
216  C           = rhoConst*SPalpha*drF/dt <-- a function of SPalpha  C           = rhoConst*SPalpha*drF(1)/dt <-- a function of SPalpha
217  C (B) SPbrinesalt provided  C (B) SPbrinesalt provided
218  C     dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt  C     dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt
219    
220  C Assuming we go with (A) here:  C Assuming we go with (B) here:
221        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
222         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
223          SPbrineSalt(i,j)= saltPlumeFlux(i,j)/SPalpha  C        SPbrineSalt(i,j,bi,bj)= saltPlumeFlux(i,j,bi,bj)/SPalpha
224       &                   *mass2rUnit*recip_drF(1)*deltaT  C     &                         *mass2rUnit*recip_drF(1)*dTtracerLev(1)
225          SPbrineSalt(i,j)=min(SPbrineSalt(i,j),SPbrineSaltmax)          SPbrineSalt(i,j,bi,bj) = SPbrineSconst
226          dMbdt(i,j) = saltPlumeFlux(i,j)/SPbrineSalt(i,j)          SPbrineSalt(i,j,bi,bj) = min( SPbrineSalt(i,j,bi,bj)
227          dVbdt(i,j) = dMbdt(i,j)*mass2rUnit       &                               ,SPbrineSaltmax )
228            dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSalt(i,j,bi,bj)
229            dVbdt(i,j)=dMbdt(i,j)*mass2rUnit
230         ENDDO         ENDDO
231        ENDDO        ENDDO
232    
# Line 237  C Distributing down: this is always from Line 234  C Distributing down: this is always from
234        DO k=Nr,1,-1        DO k=Nr,1,-1
235         DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
236          DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
237           dplumek=plumek(i,j,k+1)-plumek(i,j,k)           dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
238           dSPvolSurf2kLev(i,j,k)=dplumek*dVbdt(i,j)           dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*dVbdt(i,j)
239          ENDDO          ENDDO
240         ENDDO         ENDDO
241        ENDDO        ENDDO
# Line 246  C Distributing down: this is always from Line 243  C Distributing down: this is always from
243  C Now volume up: need to scan from bottom of SPDepth  C Now volume up: need to scan from bottom of SPDepth
244        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
245         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
246          Nlev=kLevCbot(i,j)          Nlev=SPkBottom(i,j,bi,bj)
247          IF(Nlev.GE.1) THEN          IF(Nlev.GE.1) THEN
248           DO k=Nlev,1,-1           DO k=Nlev,1,-1
249            dSPvolBelow2kLev(i,j,k) =-dSPvolkLev2Above(i,j,k+1)            dSPvolBelow2kLev(i,j,k,bi,bj)=-dSPvolkLev2Above(i,j,k+1,bi,bj)
250            dSPvolkLev2Above(i,j,k) =-( dSPvolBelow2kLev(i,j,k)            dSPvolkLev2Above(i,j,k,bi,bj)=-(dSPvolBelow2kLev(i,j,k,bi,bj)
251       &                             +  dSPvolSurf2kLev(i,j,k) )       &                                  + dSPvolSurf2kLev(i,j,k,bi,bj))
252           ENDDO           ENDDO
253          ENDIF          ENDIF
254         ENDDO         ENDDO
255        ENDDO        ENDDO
256                
257    C#ifdef ALLOW_DIAGNOSTICS
258    C      IF ( useDiagnostics ) THEN
259    C       CALL DIAGNOSTICS_FILL(
260    C     &      SPkBottom,'SPkbottm',0,1,1,bi,bj,myThid )
261    C       CALL DIAGNOSTICS_FILL(
262    C     &      SPbrineSalt,'SPbrineS',0,1,1,bi,bj,myThid )
263    C       CALL DIAGNOSTICS_FILL(
264    C     &      SPplumek,'PLUMEKB ',0,Nr,1,bi,bj,myThid )
265    C       CALL DIAGNOSTICS_FILL(
266    C     &      dSPvolSurf2kLev,'SPVs2k  ',0,Nr,1,bi,bj,myThid )
267    C       CALL DIAGNOSTICS_FILL(
268    C     &      dSPvolBelow2kLev,'SPVp2k  ',0,Nr,1,bi,bj,myThid )
269    C       CALL DIAGNOSTICS_FILL(
270    C     &      dSPvolkLev2Above,'SPVk2m  ',0,Nr,1,bi,bj,myThid )
271    C      ENDIF
272    C#endif /* ALLOW_DIAGNOSTICS */
273    
274    #endif /* SALT_PLUME_VOLUME */
275  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
276    
277        RETURN        RETURN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22