tutorial_global_oce_latlon/code/ptracers_forcing_surf.F

C $Header: /u/gcmpack/MITgcm/pkg/ptracers/ptracers_forcing_surf.F,v 1.9 2007/11/05 18:48:04 jmc Exp $
C $Name:  $

#include "PTRACERS_OPTIONS.h"

CBOP
C !ROUTINE: PTRACERS_FORCING_SURF

C !INTERFACE: ==========================================================
      SUBROUTINE PTRACERS_FORCING_SURF(
     I                            relaxForcingS,
     I                            bi, bj, iMin, iMax, jMin, jMax,
     I                            myTime,myIter,myThid )

C !DESCRIPTION:
C     Precomputes surface forcing term for pkg/ptracers.
C     Precomputation is needed because of non-local KPP transport term,
C     routine KPP_TRANSPORT_PTR.

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#include "PTRACERS_FIELDS.h"

C !INPUT PARAMETERS: ===================================================
C  relaxForcingS        :: Salt forcing due to surface relaxation
C  bi,bj                :: tile indices
C  myTime               :: model time
C  myIter               :: time-step number
C  myThid               :: thread number
      _RL relaxForcingS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER bi, bj, iMin, iMax, jMin, jMax
      _RL myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_PTRACERS

C !LOCAL VARIABLES: ====================================================
C  i,j                  :: loop indices
C  iTrc                 :: tracer index
C  ks                   :: surface level index
      INTEGER i, j
      INTEGER iTrc, ks
      _RL add2EmP(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL epsil, cutoff, tmpVar
CEOP

      IF ( usingPCoords ) THEN
        ks = Nr
      ELSE
        ks = 1
      ENDIF

C Example of how to add forcing at the surface
      DO iTrc=1,PTRACERS_numInUse
          DO j = jMin, jMax
           DO i = iMin, iMax
             surfaceForcingPTr(i,j,bi,bj,iTrc) =
     &        + 1. _d 0 / (10. _d 0 * 86400. _d 0)
     &                  * ( 0. _d 0 - pTracer(i,j,ks,bi,bj,iTrc) )
     &                  * drF(ks) * _hFacC(i,j,ks,bi,bj)
           ENDDO
          ENDDO
      ENDDO

C--   Option to convert Salt-relaxation into additional EmP contribution
      IF ( PTRACERS_addSrelax2EmP ) THEN
C-    here we assume that salt_EvPrRn = 0
C     set cutoff value to prevent too large additional EmP:
C       current limit is set to 0.1 CFL
        epsil = 1. _d -10
        cutoff = 0.1 _d 0 *drF(ks)/PTRACERS_dTLev(ks)
        IF ( ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
     &         .AND. useRealFreshWaterFlux )
     &     .OR.convertFW2Salt .EQ. -1. ) THEN
         DO j = jMin, jMax
          DO i = iMin, iMax
            tmpVar = MAX( salt(i,j,ks,bi,bj), epsil )
            add2EmP(i,j) = relaxForcingS(i,j)/tmpVar
            add2EmP(i,j) = rUnit2mass
     &                  *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
          ENDDO
         ENDDO
        ELSE
         DO j = jMin, jMax
          DO i = iMin, iMax
            add2EmP(i,j) = relaxForcingS(i,j)/convertFW2Salt
            add2EmP(i,j) = rUnit2mass
     &                  *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
          ENDDO
         ENDDO
        ENDIF
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
         CALL DIAGNOSTICS_FILL(add2EmP,'Add2EmP ',0,1,2,bi,bj,myThid)
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */
      ELSE
        DO j = jMin, jMax
          DO i = iMin, iMax
            add2EmP(i,j) = 0. _d 0
          ENDDO
        ENDDO
      ENDIF
C-- end of "addEmP" setting

#ifdef EXACT_CONSERV
      IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
     &     .AND. useRealFreshWaterFlux ) THEN

       DO iTrc=1,PTRACERS_numInUse

c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
c   the water column height ; temp., salt, (tracer) flux associated
c   with this input/output of water is added here to the surface tendency.
c
         IF (PTRACERS_EvPrRn(iTrc).NE.UNSET_RL) THEN
          DO j = jMin, jMax
           DO i = iMin, iMax
             surfaceForcingPTr(i,j,bi,bj,iTrc) =
     &          surfaceForcingPTr(i,j,bi,bj,iTrc)
     &        + ( PmEpR(i,j,bi,bj) - add2EmP(i,j) )
     &          *( PTRACERS_EvPrRn(iTrc) - pTracer(i,j,ks,bi,bj,iTrc) )
     &          *mass2rUnit
           ENDDO
          ENDDO
         ENDIF

       ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      ELSE
#else /* EXACT_CONSERV */
      IF (.TRUE.) THEN
#endif /* EXACT_CONSERV */

C--   EmPmR does not really affect the water column height (for tracer budget)
C     and is converted to a salt tendency.

       IF (convertFW2Salt .EQ. -1.) THEN
C-    use local surface tracer field to calculate forcing term:

        DO iTrc=1,PTRACERS_numInUse

         IF (PTRACERS_EvPrRn(iTrc).NE.UNSET_RL) THEN
C        account for Rain/Evap tracer content (PTRACERS_EvPrRn) using
C        local surface tracer
          DO j = jMin, jMax
           DO i = iMin, iMax
            surfaceForcingPTr(i,j,bi,bj,iTrc) =
     &          surfaceForcingPTr(i,j,bi,bj,iTrc)
     &        + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
     &          *( pTracer(i,j,ks,bi,bj,iTrc) - PTRACERS_EvPrRn(iTrc) )
     &          *mass2rUnit
           ENDDO
          ENDDO
         ENDIF

        ENDDO

       ELSE
C-    use uniform tracer value to calculate forcing term:

        DO iTrc=1,PTRACERS_numInUse

         IF (PTRACERS_EvPrRn(iTrc).NE.UNSET_RL) THEN
C     account for Rain/Evap tracer content (PTRACERS_EvPrRn) assuming uniform
C     surface tracer (=PTRACERS_ref)
          DO j = jMin, jMax
           DO i = iMin, iMax
            surfaceForcingPTr(i,j,bi,bj,iTrc) =
     &          surfaceForcingPTr(i,j,bi,bj,iTrc)
     &        + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
     &            *( PTRACERS_ref(ks,iTrc) - PTRACERS_EvPrRn(iTrc) )
     &            *mass2rUnit
           ENDDO
          ENDDO
         ENDIF

        ENDDO

C-    end local-surface-tracer / uniform-value distinction
       ENDIF

      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#endif /* ALLOW_PTRACERS */

      RETURN
      END
1	jmc	1.6	C $Header: /u/gcmpack/MITgcm/pkg/ptracers/ptracers_forcing_surf.F,v 1.9 2007/11/05 18:48:04 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "PTRACERS_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: PTRACERS_FORCING_SURF
8
9			C !INTERFACE: ==========================================================
10			SUBROUTINE PTRACERS_FORCING_SURF(
11	jmc	1.6	I relaxForcingS,
12	jmc	1.1	I bi, bj, iMin, iMax, jMin, jMax,
13			I myTime,myIter,myThid )
14
15			C !DESCRIPTION:
16			C Precomputes surface forcing term for pkg/ptracers.
17			C Precomputation is needed because of non-local KPP transport term,
18			C routine KPP_TRANSPORT_PTR.
19
20			C !USES: ===============================================================
21			IMPLICIT NONE
22			#include "SIZE.h"
23			#include "EEPARAMS.h"
24			#include "PARAMS.h"
25			#include "GRID.h"
26			#include "SURFACE.h"
27	jmc	1.6	#include "DYNVARS.h"
28	jmc	1.1	#include "FFIELDS.h"
29			#include "PTRACERS_SIZE.h"
30	jmc	1.5	#include "PTRACERS_PARAMS.h"
31			#include "PTRACERS_FIELDS.h"
32	jmc	1.1
33			C !INPUT PARAMETERS: ===================================================
34	jmc	1.6	C relaxForcingS :: Salt forcing due to surface relaxation
35	jmc	1.1	C bi,bj :: tile indices
36			C myTime :: model time
37			C myIter :: time-step number
38			C myThid :: thread number
39	jmc	1.6	_RL relaxForcingS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40	jmc	1.1	INTEGER bi, bj, iMin, iMax, jMin, jMax
41			_RL myTime
42			INTEGER myIter
43			INTEGER myThid
44
45			#ifdef ALLOW_PTRACERS
46
47			C !LOCAL VARIABLES: ====================================================
48			C i,j :: loop indices
49			C iTrc :: tracer index
50			C ks :: surface level index
51			INTEGER i, j
52			INTEGER iTrc, ks
53	jmc	1.6	_RL add2EmP(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54			_RL epsil, cutoff, tmpVar
55	jmc	1.1	CEOP
56
57			IF ( usingPCoords ) THEN
58			ks = Nr
59			ELSE
60			ks = 1
61			ENDIF
62
63			C Example of how to add forcing at the surface
64			DO iTrc=1,PTRACERS_numInUse
65			DO j = jMin, jMax
66			DO i = iMin, iMax
67	jmc	1.5	surfaceForcingPTr(i,j,bi,bj,iTrc) =
68	jmc	1.3	& + 1. _d 0 / (10. _d 0 * 86400. _d 0)
69			& * ( 0. _d 0 - pTracer(i,j,ks,bi,bj,iTrc) )
70	dfer	1.2	& * drF(ks) * _hFacC(i,j,ks,bi,bj)
71	jmc	1.1	ENDDO
72			ENDDO
73			ENDDO
74
75	jmc	1.6	C-- Option to convert Salt-relaxation into additional EmP contribution
76			IF ( PTRACERS_addSrelax2EmP ) THEN
77			C- here we assume that salt_EvPrRn = 0
78			C set cutoff value to prevent too large additional EmP:
79			C current limit is set to 0.1 CFL
80			epsil = 1. _d -10
81			cutoff = 0.1 _d 0 *drF(ks)/PTRACERS_dTLev(ks)
82			IF ( ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
83			& .AND. useRealFreshWaterFlux )
84			& .OR.convertFW2Salt .EQ. -1. ) THEN
85			DO j = jMin, jMax
86			DO i = iMin, iMax
87			tmpVar = MAX( salt(i,j,ks,bi,bj), epsil )
88			add2EmP(i,j) = relaxForcingS(i,j)/tmpVar
89			add2EmP(i,j) = rUnit2mass
90			& *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
91			ENDDO
92			ENDDO
93			ELSE
94			DO j = jMin, jMax
95			DO i = iMin, iMax
96			add2EmP(i,j) = relaxForcingS(i,j)/convertFW2Salt
97			add2EmP(i,j) = rUnit2mass
98			& *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
99			ENDDO
100			ENDDO
101			ENDIF
102			#ifdef ALLOW_DIAGNOSTICS
103			IF ( useDiagnostics ) THEN
104			CALL DIAGNOSTICS_FILL(add2EmP,'Add2EmP ',0,1,2,bi,bj,myThid)
105			ENDIF
106			#endif /* ALLOW_DIAGNOSTICS */
107			ELSE
108			DO j = jMin, jMax
109			DO i = iMin, iMax
110			add2EmP(i,j) = 0. _d 0
111			ENDDO
112			ENDDO
113			ENDIF
114			C-- end of "addEmP" setting
115
116	jmc	1.1	#ifdef EXACT_CONSERV
117			IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
118			& .AND. useRealFreshWaterFlux ) THEN
119
120			DO iTrc=1,PTRACERS_numInUse
121
122			c- NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
123			c the water column height ; temp., salt, (tracer) flux associated
124			c with this input/output of water is added here to the surface tendency.
125			c
126	jmc	1.5	IF (PTRACERS_EvPrRn(iTrc).NE.UNSET_RL) THEN
127	jmc	1.1	DO j = jMin, jMax
128			DO i = iMin, iMax
129	jmc	1.5	surfaceForcingPTr(i,j,bi,bj,iTrc) =
130			& surfaceForcingPTr(i,j,bi,bj,iTrc)
131	jmc	1.6	& + ( PmEpR(i,j,bi,bj) - add2EmP(i,j) )
132	jmc	1.1	& *( PTRACERS_EvPrRn(iTrc) - pTracer(i,j,ks,bi,bj,iTrc) )
133	jmc	1.3	& *mass2rUnit
134	jmc	1.1	ENDDO
135			ENDDO
136			ENDIF
137
138			ENDDO
139	dfer	1.4
140			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
141			ELSE
142			#else /* EXACT_CONSERV */
143			IF (.TRUE.) THEN
144			#endif /* EXACT_CONSERV */
145
146			C-- EmPmR does not really affect the water column height (for tracer budget)
147			C and is converted to a salt tendency.
148
149			IF (convertFW2Salt .EQ. -1.) THEN
150			C- use local surface tracer field to calculate forcing term:
151
152			DO iTrc=1,PTRACERS_numInUse
153
154	jmc	1.5	IF (PTRACERS_EvPrRn(iTrc).NE.UNSET_RL) THEN
155			C account for Rain/Evap tracer content (PTRACERS_EvPrRn) using
156	dfer	1.4	C local surface tracer
157			DO j = jMin, jMax
158			DO i = iMin, iMax
159	jmc	1.5	surfaceForcingPTr(i,j,bi,bj,iTrc) =
160			& surfaceForcingPTr(i,j,bi,bj,iTrc)
161	jmc	1.6	& + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
162	jmc	1.5	& *( pTracer(i,j,ks,bi,bj,iTrc) - PTRACERS_EvPrRn(iTrc) )
163	dfer	1.4	& *mass2rUnit
164			ENDDO
165			ENDDO
166			ENDIF
167
168			ENDDO
169
170			ELSE
171			C- use uniform tracer value to calculate forcing term:
172
173			DO iTrc=1,PTRACERS_numInUse
174
175	jmc	1.5	IF (PTRACERS_EvPrRn(iTrc).NE.UNSET_RL) THEN
176			C account for Rain/Evap tracer content (PTRACERS_EvPrRn) assuming uniform
177	dfer	1.4	C surface tracer (=PTRACERS_ref)
178			DO j = jMin, jMax
179			DO i = iMin, iMax
180	jmc	1.5	surfaceForcingPTr(i,j,bi,bj,iTrc) =
181			& surfaceForcingPTr(i,j,bi,bj,iTrc)
182	jmc	1.6	& + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
183	jmc	1.5	& *( PTRACERS_ref(ks,iTrc) - PTRACERS_EvPrRn(iTrc) )
184	dfer	1.4	& *mass2rUnit
185			ENDDO
186			ENDDO
187			ENDIF
188
189			ENDDO
190
191			C- end local-surface-tracer / uniform-value distinction
192			ENDIF
193
194	jmc	1.1	ENDIF
195	dfer	1.4
196			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
197	jmc	1.1
198			#endif /* ALLOW_PTRACERS */
199
200			RETURN
201			END