/[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.8 - (hide annotations) (download)
Mon Mar 5 14:52:16 2012 UTC (13 years, 4 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint64a, checkpoint64b, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63k, checkpoint64
Changes since 1.7: +2 -2 lines
- seaice_growth.F : fix pathological case treatment for SEAICE_MULTICATEGORY.
  In case ice just disappeared, reset TICES to celsius2K as done already for TICE.
- seaice_tracer_phys.F : fix typo in diagnostic name ('Fx' not 'FX').
- no change to tested results.

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

  ViewVC Help
Powered by ViewVC 1.1.22