1 |
gforget |
1.1 |
C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_tracer_growth.F,v 1.1 2011/06/02 16:21:08 gforget Exp $ |
2 |
|
|
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 |
|
|
_RL SItrFromSnow (1:sNx,1:sNy) |
41 |
|
|
_RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1 |
42 |
|
|
|
43 |
|
|
#ifdef ALLOW_SITRACER_DIAG |
44 |
|
|
_RL DIAGarray (1:sNx,1:sNy,Nr) |
45 |
|
|
#endif |
46 |
|
|
|
47 |
|
|
cgf for now I do not fully account for ocean-ice fluxes of tracer |
48 |
|
|
cgf -> I just prescribe it consistent with age tracer |
49 |
|
|
cgf eventually I will need to handle them as function params |
50 |
|
|
|
51 |
|
|
ks=1 |
52 |
|
|
|
53 |
|
|
DO bj=myByLo(myThid),myByHi(myThid) |
54 |
|
|
DO bi=myBxLo(myThid),myBxHi(myThid) |
55 |
|
|
DO iTr=1,SItrMaxNum |
56 |
|
|
|
57 |
|
|
c 0) set ice-ocean and ice-snow exchange values |
58 |
|
|
c ============================================= |
59 |
|
|
DO J=1,sNy |
60 |
|
|
DO I=1,sNx |
61 |
|
|
SItrFromOcean(i,j)=0. _d 0 |
62 |
|
|
SItrFromSnow(i,j)=0. _d 0 |
63 |
|
|
ENDDO |
64 |
|
|
ENDDO |
65 |
|
|
if (SItrName(iTr).EQ.'age') then |
66 |
|
|
c age tracer: not passed from ocean or snow |
67 |
|
|
DO J=1,sNy |
68 |
|
|
DO I=1,sNx |
69 |
|
|
SItrFromOcean(i,j)=0. _d 0 |
70 |
|
|
ENDDO |
71 |
|
|
ENDDO |
72 |
|
|
elseif (SItrName(iTr).EQ.'salinity') then |
73 |
|
|
c salinity tracer: |
74 |
|
|
DO J=1,sNy |
75 |
|
|
DO I=1,sNx |
76 |
|
|
SItrFromOcean(i,j)=SIsal0 |
77 |
|
|
#ifdef SEAICE_VARIABLE_SALINITY |
78 |
|
|
if (SIsalFRAC.GT.0.) |
79 |
|
|
& SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj) |
80 |
|
|
#endif |
81 |
|
|
ENDDO |
82 |
|
|
ENDDO |
83 |
|
|
elseif (SItrName(iTr).EQ.'one') then |
84 |
|
|
c "ice concentration" tracer that should remain .EQ.1. |
85 |
|
|
DO J=1,sNy |
86 |
|
|
DO I=1,sNx |
87 |
|
|
SItrFromOcean(i,j)=1. _d 0 |
88 |
|
|
SItrFromSnow(i,j)=1. _d 0 |
89 |
|
|
ENDDO |
90 |
|
|
ENDDO |
91 |
|
|
endif |
92 |
|
|
c 1) seaice thermodynamics processes |
93 |
|
|
c ================================== |
94 |
|
|
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 |
|
|
& +SItrFromSnow(i,j) *(1. _d 0 - growFact) |
130 |
|
|
if ((SItrName(iTr).EQ.'age').OR.(SItrName(iTr).EQ.'one')) then |
131 |
|
|
c Including this term in the bucket only makes sense for diagnostics purposes |
132 |
|
|
c and should not be done for tracers that are exchanged with the ocean (we would |
133 |
|
|
c need another bucket for ice-snow exchange, and snow tracers to start with). |
134 |
|
|
SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr) |
135 |
|
|
& -HEFFpost*SItrFromSnow(i,j)*(1. _d 0 - growFact) |
136 |
|
|
endif |
137 |
|
|
#ifdef ALLOW_SITRACER_DIAG |
138 |
|
|
diagArray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr) |
139 |
|
|
& +SItrBucket(i,j,bi,bj,iTr)-diagArray(I,J,5+(iTr-1)*5) |
140 |
|
|
#endif |
141 |
|
|
ENDDO |
142 |
|
|
ENDDO |
143 |
|
|
c 2) very ice tracer processes |
144 |
|
|
c ============================ |
145 |
|
|
if (SItrName(iTr).EQ.'age') then |
146 |
|
|
c age tracer: grow old as time passes by |
147 |
|
|
DO J=1,sNy |
148 |
|
|
DO I=1,sNx |
149 |
|
|
if (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0) then |
150 |
|
|
SItracer(i,j,bi,bj,iTr)= |
151 |
|
|
& SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm |
152 |
|
|
else |
153 |
|
|
SItracer(i,j,bi,bj,iTr)=0. _d 0 |
154 |
|
|
endif |
155 |
|
|
ENDDO |
156 |
|
|
ENDDO |
157 |
|
|
elseif (SItrName(iTr).EQ.'salinity') then |
158 |
|
|
c salinity tracer: no specific process |
159 |
|
|
elseif (SItrName(iTr).EQ.'one') then |
160 |
|
|
c "ice concentration" tracer: no specific process |
161 |
|
|
endif |
162 |
|
|
c 3) ice-ocean tracer exchange |
163 |
|
|
c ============================= |
164 |
|
|
if (SItrName(iTr).EQ.'age') then |
165 |
|
|
c age tracer: not passed to ocean |
166 |
|
|
elseif (SItrName(iTr).EQ.'salinity') then |
167 |
|
|
c salinity tracer: salt flux |
168 |
|
|
c DO J=1,sNy |
169 |
|
|
c DO I=1,sNx |
170 |
|
|
c saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr) |
171 |
|
|
c & *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce |
172 |
|
|
c note: at this point of the time step, that is the correct sign |
173 |
|
|
c saltPlumeFlux(I,J,bi,bj) = ... |
174 |
|
|
c ENDDO |
175 |
|
|
c ENDDO |
176 |
|
|
elseif (SItrName(iTr).EQ.'one') then |
177 |
|
|
c "ice concentration" tracer: not passed to ocean |
178 |
|
|
endif |
179 |
|
|
DO J=1,sNy |
180 |
|
|
DO I=1,sNx |
181 |
|
|
#ifdef ALLOW_SITRACER_DIAG |
182 |
|
|
diagArray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr) |
183 |
|
|
& *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce |
184 |
|
|
#endif |
185 |
|
|
c empty bucket |
186 |
|
|
SItrBucket(i,j,bi,bj,iTr)=0. _d 0 |
187 |
|
|
ENDDO |
188 |
|
|
ENDDO |
189 |
|
|
c 4) diagnostics |
190 |
|
|
c ============== |
191 |
|
|
#ifdef ALLOW_SITRACER_DIAG |
192 |
|
|
DO J=1,sNy |
193 |
|
|
DO I=1,sNx |
194 |
|
|
HEFFpost=SItrHEFF(i,j,bi,bj,5) |
195 |
|
|
DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr) |
196 |
|
|
DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*HEFFpost |
197 |
|
|
c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2) |
198 |
|
|
if (SItrName(iTr).EQ.'age') then |
199 |
|
|
DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2) |
200 |
|
|
elseif (SItrName(iTr).EQ.'salinity') then |
201 |
|
|
DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce |
202 |
|
|
elseif (SItrName(iTr).EQ.'one') then |
203 |
|
|
DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost |
204 |
|
|
endif |
205 |
|
|
c diagArray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys |
206 |
|
|
c diagArray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer) |
207 |
|
|
ENDDO |
208 |
|
|
ENDDO |
209 |
|
|
#endif |
210 |
|
|
ENDDO |
211 |
|
|
#ifdef ALLOW_SITRACER_DIAG |
212 |
|
|
CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1 ',0,Nr,3,bi,bj,myThid) |
213 |
|
|
#endif |
214 |
|
|
ENDDO |
215 |
|
|
ENDDO |
216 |
|
|
|
217 |
|
|
#endif /* ALLOW_SITRACER */ |
218 |
|
|
|
219 |
|
|
RETURN |
220 |
|
|
END |