1 |
gforget |
1.8 |
C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_tracer_phys.F,v 1.7 2012/02/16 03:32:39 jmc Exp $ |
2 |
gforget |
1.1 |
C $Name: $ |
3 |
|
|
|
4 |
|
|
#include "SEAICE_OPTIONS.h" |
5 |
|
|
|
6 |
|
|
CStartOfInterface |
7 |
|
|
SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid ) |
8 |
jmc |
1.7 |
C *=======================================================* |
9 |
|
|
C | SUBROUTINE seaice_tracer_phys |
10 |
|
|
C | o Time step SItr/SItrEFF as a result of |
11 |
|
|
C | seaice thermodynamics and specific tracer physics |
12 |
|
|
C *=======================================================* |
13 |
gforget |
1.1 |
IMPLICIT NONE |
14 |
|
|
|
15 |
|
|
C === Global variables === |
16 |
|
|
#include "SIZE.h" |
17 |
|
|
#include "EEPARAMS.h" |
18 |
|
|
#include "FFIELDS.h" |
19 |
|
|
#include "DYNVARS.h" |
20 |
|
|
#include "SEAICE_SIZE.h" |
21 |
|
|
#include "SEAICE.h" |
22 |
|
|
#include "SEAICE_PARAMS.h" |
23 |
|
|
#include "SEAICE_TRACER.h" |
24 |
gforget |
1.6 |
#ifdef ALLOW_SALT_PLUME |
25 |
|
|
# include "SALT_PLUME.h" |
26 |
|
|
#endif |
27 |
gforget |
1.1 |
|
28 |
|
|
C === Routine arguments === |
29 |
|
|
C INPUT: |
30 |
|
|
C myTime :: Simulation time |
31 |
|
|
C myIter :: Simulation timestep number |
32 |
|
|
C myThid :: Thread no. that called this routine. |
33 |
|
|
C OUTPUT: |
34 |
|
|
_RL myTime |
35 |
|
|
INTEGER myIter, myThid |
36 |
|
|
CEndOfInterface |
37 |
|
|
|
38 |
|
|
C === Local variables === |
39 |
|
|
#ifdef ALLOW_SITRACER |
40 |
|
|
|
41 |
|
|
INTEGER iTr, jTh, I, J, bi, bj, ks |
42 |
|
|
_RL SItrFromOcean (1:sNx,1:sNy) |
43 |
gforget |
1.2 |
_RL SItrFromFlood (1:sNx,1:sNy) |
44 |
gforget |
1.1 |
_RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1 |
45 |
gforget |
1.2 |
_RL SItrExpand (1:sNx,1:sNy) |
46 |
|
|
_RL AREAprev, AREApost, expandFact |
47 |
gforget |
1.3 |
CHARACTER*8 diagName |
48 |
gforget |
1.1 |
|
49 |
gforget |
1.3 |
#ifdef ALLOW_SITRACER_DEBUG_DIAG |
50 |
gforget |
1.1 |
_RL DIAGarray (1:sNx,1:sNy,Nr) |
51 |
|
|
#endif |
52 |
|
|
|
53 |
|
|
cgf for now I do not fully account for ocean-ice fluxes of tracer |
54 |
|
|
cgf -> I just prescribe it consistent with age tracer |
55 |
|
|
cgf eventually I will need to handle them as function params |
56 |
|
|
|
57 |
|
|
ks=1 |
58 |
|
|
|
59 |
|
|
DO bj=myByLo(myThid),myByHi(myThid) |
60 |
|
|
DO bi=myBxLo(myThid),myBxHi(myThid) |
61 |
gforget |
1.6 |
DO iTr=1,SItrNumInUse |
62 |
gforget |
1.1 |
|
63 |
jmc |
1.4 |
c 0) set ice-ocean and ice-snow exchange values |
64 |
gforget |
1.1 |
c ============================================= |
65 |
|
|
DO J=1,sNy |
66 |
|
|
DO I=1,sNx |
67 |
gforget |
1.6 |
SItrFromOcean(i,j)=SItrFromOcean0(iTr) |
68 |
|
|
SItrFromFlood(i,j)=SItrFromFlood0(iTr) |
69 |
|
|
SItrExpand(i,j)=SItrExpand0(iTr) |
70 |
gforget |
1.1 |
ENDDO |
71 |
|
|
ENDDO |
72 |
jmc |
1.4 |
c salinity tracer: |
73 |
jmc |
1.7 |
if ( (SItrName(iTr).EQ.'salinity').AND. |
74 |
gforget |
1.6 |
& (SItrFromOceanFrac(iTr).GT.ZERO) ) then |
75 |
gforget |
1.1 |
DO J=1,sNy |
76 |
|
|
DO I=1,sNx |
77 |
gforget |
1.6 |
SItrFromOcean(i,j)=SItrFromOceanFrac(iTr)*salt(I,j,ks,bi,bj) |
78 |
|
|
SItrFromFlood(i,j)=SItrFromFloodFrac(iTr)*salt(I,j,ks,bi,bj) |
79 |
gforget |
1.1 |
ENDDO |
80 |
|
|
ENDDO |
81 |
|
|
endif |
82 |
|
|
c 1) seaice thermodynamics processes |
83 |
|
|
c ================================== |
84 |
gforget |
1.2 |
if (SItrMate(iTr).EQ.'HEFF') then |
85 |
gforget |
1.1 |
DO J=1,sNy |
86 |
|
|
DO I=1,sNx |
87 |
|
|
HEFFprev=SItrHEFF(i,j,bi,bj,1) |
88 |
gforget |
1.3 |
#ifdef ALLOW_SITRACER_DEBUG_DIAG |
89 |
jmc |
1.4 |
DIAGarray(I,J,5+(iTr-1)*5) = |
90 |
gforget |
1.1 |
& HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr) |
91 |
|
|
#endif |
92 |
|
|
c apply the sequence of thermodynamics increments to actual traceur |
93 |
|
|
c (see seaice_growth.F) |
94 |
|
|
c (jTh=1 tendency due to ice-ocean interaction) |
95 |
|
|
c (jTh=2 tendency due to the atmosphere, over ice covered part) |
96 |
|
|
c (jTh=3 tendency due to the atmosphere, over open water part) |
97 |
|
|
c (jTh=4 tendency due to flooding) |
98 |
|
|
DO jTh=1,3 |
99 |
|
|
HEFFprev=SItrHEFF(i,j,bi,bj,jTh) |
100 |
|
|
HEFFpost=SItrHEFF(i,j,bi,bj,jTh+1) |
101 |
|
|
c compute ratio in [0. 1.] range for either growth or melt |
102 |
|
|
growFact=1. _d 0 |
103 |
|
|
meltPart=0. _d 0 |
104 |
|
|
if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost |
105 |
|
|
if (HEFFpost.LT.HEFFprev) meltPart=HEFFprev-HEFFpost |
106 |
|
|
c update SItr accordingly |
107 |
|
|
SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact |
108 |
|
|
& +SItrFromOcean(i,j)*(1. _d 0 - growFact) |
109 |
|
|
SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr) |
110 |
|
|
& -HEFFpost*SItrFromOcean(i,j)*(1. _d 0 - growFact) |
111 |
|
|
SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr) |
112 |
|
|
& +meltPart*SItracer(i,j,bi,bj,iTr) |
113 |
|
|
ENDDO |
114 |
|
|
c apply flooding term |
115 |
|
|
growFact=1. _d 0 |
116 |
|
|
HEFFprev=SItrHEFF(i,j,bi,bj,4) |
117 |
|
|
HEFFpost=SItrHEFF(i,j,bi,bj,5) |
118 |
|
|
if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost |
119 |
|
|
SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact |
120 |
gforget |
1.2 |
& +SItrFromFlood(i,j) *(1. _d 0 - growFact) |
121 |
|
|
c rk: flooding can only imply an ocean-ice tracer exchange, as long |
122 |
|
|
c as we dont have snow tracers, so it goes through SItrBucket. |
123 |
gforget |
1.1 |
SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr) |
124 |
gforget |
1.2 |
& -HEFFpost*SItrFromFlood(i,j)*(1. _d 0 - growFact) |
125 |
gforget |
1.3 |
#ifdef ALLOW_SITRACER_DEBUG_DIAG |
126 |
jmc |
1.4 |
DIAGarray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr) |
127 |
gforget |
1.3 |
& +SItrBucket(i,j,bi,bj,iTr)-DIAGarray(I,J,5+(iTr-1)*5) |
128 |
gforget |
1.1 |
#endif |
129 |
|
|
ENDDO |
130 |
|
|
ENDDO |
131 |
gforget |
1.2 |
c TAF? if (SItrMate(iTr).EQ.'AREA') then |
132 |
|
|
else |
133 |
|
|
c 1) or seaice cover expansion |
134 |
|
|
c ============================ |
135 |
jmc |
1.4 |
c this is much simpler than for ice volume/mass tracers, because |
136 |
|
|
c properties of the ice surface are not be conserved across the |
137 |
gforget |
1.2 |
c ocean-ice system, the contraction/expansion terms are all |
138 |
jmc |
1.4 |
c simultaneous (which is sane), and the only generic effect |
139 |
gforget |
1.2 |
c is due to expansion (new cover). |
140 |
|
|
DO J=1,sNy |
141 |
|
|
DO I=1,sNx |
142 |
|
|
c apply expansion |
143 |
|
|
AREAprev=SItrAREA(i,j,bi,bj,2) |
144 |
|
|
AREApost=SItrAREA(i,j,bi,bj,3) |
145 |
|
|
c compute ratio in [0. 1.] range for expansion/contraction |
146 |
|
|
expandFact=1. _d 0 |
147 |
|
|
if (AREApost.GT.AREAprev) expandFact=AREAprev/AREApost |
148 |
|
|
c update SItr accordingly |
149 |
|
|
SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*expandFact |
150 |
|
|
& +SItrExpand(i,j)*(1. _d 0 - expandFact) |
151 |
|
|
ENDDO |
152 |
|
|
ENDDO |
153 |
|
|
endif |
154 |
gforget |
1.1 |
c 2) very ice tracer processes |
155 |
|
|
c ============================ |
156 |
|
|
if (SItrName(iTr).EQ.'age') then |
157 |
|
|
c age tracer: grow old as time passes by |
158 |
|
|
DO J=1,sNy |
159 |
|
|
DO I=1,sNx |
160 |
gforget |
1.2 |
if (( (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0).AND.(SItrMate(iTr) |
161 |
|
|
& .EQ.'HEFF') ).OR.( (SItrAREA(i,j,bi,bj,3).GT.0. _d 0).AND. |
162 |
|
|
& (SItrMate(iTr).EQ.'AREA') )) then |
163 |
gforget |
1.1 |
SItracer(i,j,bi,bj,iTr)= |
164 |
|
|
& SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm |
165 |
|
|
else |
166 |
|
|
SItracer(i,j,bi,bj,iTr)=0. _d 0 |
167 |
|
|
endif |
168 |
|
|
ENDDO |
169 |
|
|
ENDDO |
170 |
|
|
elseif (SItrName(iTr).EQ.'salinity') then |
171 |
|
|
c salinity tracer: no specific process |
172 |
|
|
elseif (SItrName(iTr).EQ.'one') then |
173 |
|
|
c "ice concentration" tracer: no specific process |
174 |
gforget |
1.2 |
elseif (SItrName(iTr).EQ.'ridge') then |
175 |
|
|
c simple, made up, ice surface roughness index prototype |
176 |
|
|
#ifndef SEAICE_GROWTH_LEGACY |
177 |
|
|
DO J=1,sNy |
178 |
|
|
DO I=1,sNx |
179 |
|
|
c ridging increases roughness |
180 |
|
|
SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)+ |
181 |
|
|
& MAX(0. _d 0, SItrAREA(i,j,bi,bj,1)-SItrAREA(i,j,bi,bj,2)) |
182 |
|
|
c ice melt reduces ridges/roughness |
183 |
|
|
HEFFprev=SItrHEFF(i,j,bi,bj,1) |
184 |
|
|
HEFFpost=SItrHEFF(i,j,bi,bj,4) |
185 |
|
|
tmpscal1=1. _d 0 |
186 |
|
|
if (HEFFprev.GT.HEFFpost) tmpscal1=HEFFpost/HEFFprev |
187 |
|
|
SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*tmpscal1 |
188 |
|
|
ENDDO |
189 |
|
|
ENDDO |
190 |
|
|
#endif |
191 |
gforget |
1.1 |
endif |
192 |
gforget |
1.2 |
c 3) ice-ocean tracer exchange/mapping to external variables |
193 |
|
|
c ========================================================== |
194 |
gforget |
1.3 |
#ifdef ALLOW_DIAGNOSTICS |
195 |
|
|
if (SItrMate(iTr).EQ.'HEFF') then |
196 |
gforget |
1.8 |
WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'Fx' |
197 |
gforget |
1.3 |
tmpscal1=-ONE/SEAICE_deltaTtherm*SEAICE_rhoIce |
198 |
jmc |
1.7 |
CALL DIAGNOSTICS_SCALE_FILL(SItrBucket(1-OLx,1-OLy,bi,bj,iTr), |
199 |
gforget |
1.3 |
& tmpscal1, 1, diagName,0,1,2,bi,bj,myThid) |
200 |
|
|
endif |
201 |
|
|
#endif |
202 |
gforget |
1.6 |
|
203 |
|
|
if ( (SItrName(iTr).EQ.'salinity').AND. |
204 |
|
|
& (SEAICE_salinityTracer) ) then |
205 |
gforget |
1.1 |
c salinity tracer: salt flux |
206 |
gforget |
1.6 |
DO J=1,sNy |
207 |
|
|
DO I=1,sNx |
208 |
|
|
saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr) |
209 |
|
|
& *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce |
210 |
gforget |
1.1 |
c note: at this point of the time step, that is the correct sign |
211 |
gforget |
1.6 |
#ifdef ALLOW_SALT_PLUME |
212 |
|
|
c should work for both constant and variable ice salinity -- to be tested |
213 |
jmc |
1.7 |
saltPlumeFlux(I,J,bi,bj) = MAX(ZERO,saltFlux(I,J,bi,bj)) |
214 |
gforget |
1.6 |
& *SPsalFRAC*(salt(I,j,ks,bi,bj)-SItrFromOcean(i,j)) |
215 |
|
|
#endif |
216 |
|
|
ENDDO |
217 |
|
|
ENDDO |
218 |
gforget |
1.1 |
endif |
219 |
jmc |
1.7 |
|
220 |
gforget |
1.1 |
DO J=1,sNy |
221 |
|
|
DO I=1,sNx |
222 |
gforget |
1.3 |
#ifdef ALLOW_SITRACER_DEBUG_DIAG |
223 |
|
|
DIAGarray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr) |
224 |
gforget |
1.1 |
& *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce |
225 |
|
|
#endif |
226 |
|
|
c empty bucket |
227 |
|
|
SItrBucket(i,j,bi,bj,iTr)=0. _d 0 |
228 |
|
|
ENDDO |
229 |
|
|
ENDDO |
230 |
jmc |
1.7 |
|
231 |
gforget |
1.2 |
c TAF? elseif (SItrMate(iTr).EQ.'AREA') then |
232 |
gforget |
1.6 |
|
233 |
gforget |
1.1 |
c 4) diagnostics |
234 |
|
|
c ============== |
235 |
gforget |
1.3 |
#ifdef ALLOW_SITRACER_DEBUG_DIAG |
236 |
gforget |
1.2 |
if (SItrMate(iTr).EQ.'HEFF') then |
237 |
gforget |
1.1 |
DO J=1,sNy |
238 |
|
|
DO I=1,sNx |
239 |
|
|
HEFFpost=SItrHEFF(i,j,bi,bj,5) |
240 |
|
|
DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr) |
241 |
|
|
DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*HEFFpost |
242 |
gforget |
1.3 |
c DIAGarray(:,:,3) is the term of comparison for DIAGarray(:,:,2) |
243 |
gforget |
1.5 |
if (SItrName(iTr).EQ.'salinity') then |
244 |
gforget |
1.1 |
DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce |
245 |
|
|
elseif (SItrName(iTr).EQ.'one') then |
246 |
|
|
DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost |
247 |
|
|
endif |
248 |
gforget |
1.3 |
c DIAGarray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys |
249 |
|
|
c DIAGarray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer) |
250 |
gforget |
1.1 |
ENDDO |
251 |
|
|
ENDDO |
252 |
gforget |
1.2 |
else |
253 |
|
|
DO J=1,sNy |
254 |
|
|
DO I=1,sNx |
255 |
|
|
AREApost=SItrAREA(i,j,bi,bj,3) |
256 |
|
|
DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr) |
257 |
|
|
DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*AREApost |
258 |
|
|
ENDDO |
259 |
|
|
ENDDO |
260 |
|
|
endif |
261 |
gforget |
1.1 |
#endif |
262 |
|
|
ENDDO |
263 |
gforget |
1.3 |
#ifdef ALLOW_SITRACER_DEBUG_DIAG |
264 |
jmc |
1.4 |
c CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1 ',0,Nr,3,bi,bj,myThid) |
265 |
gforget |
1.1 |
#endif |
266 |
|
|
ENDDO |
267 |
|
|
ENDDO |
268 |
|
|
|
269 |
|
|
#endif /* ALLOW_SITRACER */ |
270 |
|
|
|
271 |
|
|
RETURN |
272 |
|
|
END |