atnguyen/code_21Dec2012_saltplume/salt_plume_tendency_apply_s.F

C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_tendency_apply_s.F,v 1.3 2014/04/20 09:00:36 atn Exp $
C $Name:  $

#include "SALT_PLUME_OPTIONS.h"

CBOP 0
C     !ROUTINE: SALT_PLUME_TENDENCY_APPLY_S
C     !INTERFACE:
      SUBROUTINE SALT_PLUME_TENDENCY_APPLY_S(
     &     iMin, iMax, jMin, jMax,
     &     bi,bj,kLev,myTime,myThid)

C     !DESCRIPTION:
C     Add salt_plume tendency terms to S tendency.
C     Routine works for one level at a time.
C     SaltPlume is the amount of salt rejected by ice while freezing;
C     it is here redistributed to multiple vertical levels as per
C     Duffy et al. (GRL 1999).

C     !INPUT PARAMETERS:
      IMPLICIT NONE
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "SALT_PLUME.h"

C     !INPUT PARAMETERS:
      integer iMin, iMax, jMin, jMax, kLev, bi, bj, myThid
      _RL myTime
CEOP

#ifdef ALLOW_SALT_PLUME
#ifndef SALT_PLUME_VOLUME

C     !LOCAL VARIABLES:
      integer i, j
      _RL minusone
      parameter(minusone = -1.)
      _RL plumefrac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL plumetend(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef TARGET_NEC_SX
      integer imt
      parameter( imt=(sNx+2*OLx)*(sNy+2*OLy) )
      _RL plumekb2D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#else
      integer two2
      parameter(two2 = 2)
      _RL plumekb(two2), SPdepth(two2)
#endif

#ifdef TARGET_NEC_SX
C     The vector version computes plumekb2D at each grid point, but this
C     is still faster than non-vector code.
      IF ( kLev .LT. Nr ) THEN
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         plumekb2D(i,j)=abs(rF(kLev))
        ENDDO
       ENDDO
       CALL SALT_PLUME_FRAC(
     I      imt,minusone,SaltPlumeDepth(1-Olx,1-Oly,bi,bj),
     U      plumekb2D,
     I      myTime, 1, myThid )
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         plumefrac(I,J) = plumekb2D(i,j)
         plumekb2D(i,j) = abs(rF(kLev+1))
        ENDDO
       ENDDO
       CALL SALT_PLUME_FRAC(
     I      imt,minusone,SaltPlumeDepth(1-Olx,1-Oly,bi,bj),
     U      plumekb2D,
     I      myTime, 1, myThid )
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
C     Penetrating saltplume fraction:cumSP(klev+1)-cumSP(klev)
         IF ( SaltPlumeDepth(i,j,bi,bj) .GT. abs(rF(kLev)) ) THEN
          plumefrac(i,j) = ( plumekb2D(i,j)-plumefrac(i,j) )
     &                     *maskC(i,j,klev,bi,bj)
          plumetend(I,J) = saltPlumeFlux(i,j,bi,bj)*plumefrac(I,J)
          gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)+plumetend(I,J)
     &        *recip_drF(kLev)*mass2rUnit*_recip_hFacC(i,j,kLev,bi,bj)
         ELSE
          plumefrac(i,j) = 0. _d 0
          plumetend(I,J) = 0. _d 0
         ENDIF
        ENDDO
       ENDDO
      ENDIF
#else
      DO j=jMin,jMax
       DO i=iMin,iMax
C Penetrating saltplume fraction:cumulativeSP(klev+1)-cumulativeSP(klev)
        IF ( SaltPlumeDepth(i,j,bi,bj) .GT. abs(rF(kLev)) ) THEN
         plumefrac(I,J) = 0. _d 0
         plumekb(1)=abs(rF(klev))
         plumekb(2)=abs(rF(klev+1))
         SPdepth(1)=SaltPlumeDepth(i,j,bi,bj)
         SPdepth(2)=SaltPlumeDepth(i,j,bi,bj)
         CALL SALT_PLUME_FRAC(
     I                   two2,minusone,SPdepth,
     U                   plumekb,
     I                   myTime, 1, myThid )
         plumefrac(I,J) = (plumekb(2)-plumekb(1))*maskC(i,j,klev,bi,bj)
         plumetend(I,J) = saltPlumeFlux(i,j,bi,bj)*plumefrac(I,J)
         gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)+plumetend(I,J)
     &        *recip_drF(kLev)*mass2rUnit*_recip_hFacC(i,j,kLev,bi,bj)
        ELSE
         plumefrac(I,J) = 0. _d 0
         plumetend(I,J) = 0. _d 0
        ENDIF
       ENDDO
      ENDDO
#endif /* TARGET_NEC_SX */

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
       CALL DIAGNOSTICS_FILL (
     &      plumefrac,'PLUMEKB ',kLev,1,2,bi,bj,myThid )
       CALL DIAGNOSTICS_FILL (
     &      plumetend,'oceSPtnd',kLev,1,2,bi,bj,myThid )
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */
      
#endif /* SALT_PLUME_VOLUME */
#endif /* ALLOW_SALT_PLUME */
      
      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_tendency_apply_s.F,v 1.3 2014/04/20 09:00:36 atn Exp $
2	C $Name: $
3
4	#include "SALT_PLUME_OPTIONS.h"
5
6	CBOP 0
7	C !ROUTINE: SALT_PLUME_TENDENCY_APPLY_S
8	C !INTERFACE:
9	SUBROUTINE SALT_PLUME_TENDENCY_APPLY_S(
10	& iMin, iMax, jMin, jMax,
11	& bi,bj,kLev,myTime,myThid)
12
13	C !DESCRIPTION:
14	C Add salt_plume tendency terms to S tendency.
15	C Routine works for one level at a time.
16	C SaltPlume is the amount of salt rejected by ice while freezing;
17	C it is here redistributed to multiple vertical levels as per
18	C Duffy et al. (GRL 1999).
19
20	C !INPUT PARAMETERS:
21	IMPLICIT NONE
22	#include "SIZE.h"
23	#include "GRID.h"
24	#include "EEPARAMS.h"
25	#include "PARAMS.h"
26	#include "DYNVARS.h"
27	#include "SALT_PLUME.h"
28
29	C !INPUT PARAMETERS:
30	integer iMin, iMax, jMin, jMax, kLev, bi, bj, myThid
31	_RL myTime
32	CEOP
33
34	#ifdef ALLOW_SALT_PLUME
35	#ifndef SALT_PLUME_VOLUME
36
37	C !LOCAL VARIABLES:
38	integer i, j
39	_RL minusone
40	parameter(minusone = -1.)
41	_RL plumefrac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42	_RL plumetend(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43	#ifdef TARGET_NEC_SX
44	integer imt
45	parameter( imt=(sNx+2OLx)(sNy+2*OLy) )
46	_RL plumekb2D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47	#else
48	integer two2
49	parameter(two2 = 2)
50	_RL plumekb(two2), SPdepth(two2)
51	#endif
52
53	#ifdef TARGET_NEC_SX
54	C The vector version computes plumekb2D at each grid point, but this
55	C is still faster than non-vector code.
56	IF ( kLev .LT. Nr ) THEN
57	DO j=1-Oly,sNy+Oly
58	DO i=1-Olx,sNx+Olx
59	plumekb2D(i,j)=abs(rF(kLev))
60	ENDDO
61	ENDDO
62	CALL SALT_PLUME_FRAC(
63	I imt,minusone,SaltPlumeDepth(1-Olx,1-Oly,bi,bj),
64	U plumekb2D,
65	I myTime, 1, myThid )
66	DO j=1-Oly,sNy+Oly
67	DO i=1-Olx,sNx+Olx
68	plumefrac(I,J) = plumekb2D(i,j)
69	plumekb2D(i,j) = abs(rF(kLev+1))
70	ENDDO
71	ENDDO
72	CALL SALT_PLUME_FRAC(
73	I imt,minusone,SaltPlumeDepth(1-Olx,1-Oly,bi,bj),
74	U plumekb2D,
75	I myTime, 1, myThid )
76	DO j=1-Oly,sNy+Oly
77	DO i=1-Olx,sNx+Olx
78	C Penetrating saltplume fraction:cumSP(klev+1)-cumSP(klev)
79	IF ( SaltPlumeDepth(i,j,bi,bj) .GT. abs(rF(kLev)) ) THEN
80	plumefrac(i,j) = ( plumekb2D(i,j)-plumefrac(i,j) )
81	& *maskC(i,j,klev,bi,bj)
82	plumetend(I,J) = saltPlumeFlux(i,j,bi,bj)*plumefrac(I,J)
83	gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)+plumetend(I,J)
84	& recip_drF(kLev)mass2rUnit*_recip_hFacC(i,j,kLev,bi,bj)
85	ELSE
86	plumefrac(i,j) = 0. _d 0
87	plumetend(I,J) = 0. _d 0
88	ENDIF
89	ENDDO
90	ENDDO
91	ENDIF
92	#else
93	DO j=jMin,jMax
94	DO i=iMin,iMax
95	C Penetrating saltplume fraction:cumulativeSP(klev+1)-cumulativeSP(klev)
96	IF ( SaltPlumeDepth(i,j,bi,bj) .GT. abs(rF(kLev)) ) THEN
97	plumefrac(I,J) = 0. _d 0
98	plumekb(1)=abs(rF(klev))
99	plumekb(2)=abs(rF(klev+1))
100	SPdepth(1)=SaltPlumeDepth(i,j,bi,bj)
101	SPdepth(2)=SaltPlumeDepth(i,j,bi,bj)
102	CALL SALT_PLUME_FRAC(
103	I two2,minusone,SPdepth,
104	U plumekb,
105	I myTime, 1, myThid )
106	plumefrac(I,J) = (plumekb(2)-plumekb(1))*maskC(i,j,klev,bi,bj)
107	plumetend(I,J) = saltPlumeFlux(i,j,bi,bj)*plumefrac(I,J)
108	gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)+plumetend(I,J)
109	& recip_drF(kLev)mass2rUnit*_recip_hFacC(i,j,kLev,bi,bj)
110	ELSE
111	plumefrac(I,J) = 0. _d 0
112	plumetend(I,J) = 0. _d 0
113	ENDIF
114	ENDDO
115	ENDDO
116	#endif /* TARGET_NEC_SX */
117
118	#ifdef ALLOW_DIAGNOSTICS
119	IF ( useDiagnostics ) THEN
120	CALL DIAGNOSTICS_FILL (
121	& plumefrac,'PLUMEKB ',kLev,1,2,bi,bj,myThid )
122	CALL DIAGNOSTICS_FILL (
123	& plumetend,'oceSPtnd',kLev,1,2,bi,bj,myThid )
124	ENDIF
125	#endif /* ALLOW_DIAGNOSTICS */
126
127	#endif /* SALT_PLUME_VOLUME */
128	#endif /* ALLOW_SALT_PLUME */
129
130	RETURN
131	END