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