/[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.2 - (hide annotations) (download)
Thu Jun 9 19:37:01 2011 UTC (14 years, 1 month ago) by gforget
Branch: MAIN
Changes since 1.1: +74 -20 lines
- seaice_tracer_phys.F and seaice_advdiff.F
  do ice cover tracers, in addition to ice volume tracers.
- seaice_readparms.F
  add SItrMate ('HEFF' or 'AREA') in PARAMS03 to switch
    from ice volume tracer (defualt) to ice cover tracer
- seaice_growth.F
  use areaMax in AREA update (part 4), consistent with ridging step (part 2.5).
  store AREA in SItrAREA at the ridging and update steps (for use with SItracer).

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

  ViewVC Help
Powered by ViewVC 1.1.22