forked from darwinproject/darwin3
-
Notifications
You must be signed in to change notification settings - Fork 1
/
ptracers_forcing_surf.F
202 lines (173 loc) · 6.15 KB
/
ptracers_forcing_surf.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#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_START.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
c IF ( PTRACERS_StepFwd(iTrc) ) THEN
DO j = jMin, jMax
DO i = iMin, iMax
surfaceForcingPTr(i,j,bi,bj,iTrc) = 0. _d 0
& + surfaceForcingS(i,j,bi,bj)
ENDDO
ENDDO
c ENDIF
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_StepFwd(iTrc) .AND.
& 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_StepFwd(iTrc) .AND.
& 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_StepFwd(iTrc) .AND.
& 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