/[MITgcm]/MITgcm/pkg/seaice/seaice_tracer_phys.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/seaice_tracer_phys.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.5 - (hide annotations) (download)
Fri Feb 3 13:34:32 2012 UTC (13 years, 5 months ago) by gforget
Branch: MAIN
Changes since 1.4: +2 -8 lines
- removal of the old way of seaice age tracer, which is now replaced by particular cases of SITRACER.
- retired params : SEAICEadvAge, SEAICEadvSchAge, SEAICEdiffKhAge, IceAgeTrFile.
- added to SITRACER : IceAgeTrFile, check pickups, monitor, output.

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

  ViewVC Help
Powered by ViewVC 1.1.22