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

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

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

revision 1.1 by gforget, Tue Jun 7 03:58:23 2011 UTC revision 1.3 by gforget, Mon Jun 13 23:21:18 2011 UTC
# Line 37  C     === Local variables === Line 37  C     === Local variables ===
37    
38        INTEGER iTr, jTh, I, J, bi, bj, ks        INTEGER iTr, jTh, I, J, bi, bj, ks
39        _RL SItrFromOcean  (1:sNx,1:sNy)        _RL SItrFromOcean  (1:sNx,1:sNy)
40        _RL SItrFromSnow   (1:sNx,1:sNy)        _RL SItrFromFlood   (1:sNx,1:sNy)
41        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
42          _RL SItrExpand  (1:sNx,1:sNy)
43          _RL AREAprev, AREApost, expandFact
44          CHARACTER*8   diagName
45    
46  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
47        _RL DIAGarray     (1:sNx,1:sNy,Nr)        _RL DIAGarray     (1:sNx,1:sNy,Nr)
48  #endif  #endif
49    
# Line 59  c ====================================== Line 62  c ======================================
62        DO J=1,sNy        DO J=1,sNy
63         DO I=1,sNx         DO I=1,sNx
64          SItrFromOcean(i,j)=0. _d 0          SItrFromOcean(i,j)=0. _d 0
65          SItrFromSnow(i,j)=0. _d 0          SItrFromFlood(i,j)=0. _d 0
66            SItrExpand(i,j)=0. _d 0
67         ENDDO         ENDDO
68        ENDDO        ENDDO
69        if (SItrName(iTr).EQ.'age') then        if (SItrName(iTr).EQ.'age') then
70  c age tracer: not passed from ocean or snow  c age tracer: no age in ocean, or effect from ice cover changes
        DO J=1,sNy  
         DO I=1,sNx  
          SItrFromOcean(i,j)=0. _d 0  
         ENDDO  
        ENDDO  
71        elseif (SItrName(iTr).EQ.'salinity') then        elseif (SItrName(iTr).EQ.'salinity') then
72  c salinity tracer:  c salinity tracer:
73         DO J=1,sNy         DO J=1,sNy
# Line 78  c salinity tracer: Line 77  c salinity tracer:
77           if (SIsalFRAC.GT.0.)           if (SIsalFRAC.GT.0.)
78       &   SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj)       &   SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj)
79  #endif  #endif
80    c as of now, flooding implies no salt extraction from ocean
81          ENDDO          ENDDO
82         ENDDO         ENDDO
83        elseif (SItrName(iTr).EQ.'one') then        elseif (SItrName(iTr).EQ.'one') then
# Line 85  c "ice concentration" tracer that should Line 85  c "ice concentration" tracer that should
85         DO J=1,sNy         DO J=1,sNy
86          DO I=1,sNx          DO I=1,sNx
87           SItrFromOcean(i,j)=1. _d 0           SItrFromOcean(i,j)=1. _d 0
88           SItrFromSnow(i,j)=1. _d 0           SItrFromFlood(i,j)=1. _d 0
89          ENDDO          ENDDO
90         ENDDO         ENDDO
91        endif        endif
92  c 1) seaice thermodynamics processes  c 1) seaice thermodynamics processes
93  c ==================================  c ==================================
94          if (SItrMate(iTr).EQ.'HEFF') then
95        DO J=1,sNy        DO J=1,sNy
96         DO I=1,sNx         DO I=1,sNx
97          HEFFprev=SItrHEFF(i,j,bi,bj,1)          HEFFprev=SItrHEFF(i,j,bi,bj,1)
98  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
99          diagArray(I,J,5+(iTr-1)*5) =          DIAGarray(I,J,5+(iTr-1)*5) =
100       &    HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr)       &    HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr)
101  #endif  #endif
102  c apply the sequence of thermodynamics increments to actual traceur  c apply the sequence of thermodynamics increments to actual traceur
# Line 126  c apply flooding term Line 127  c apply flooding term
127          HEFFpost=SItrHEFF(i,j,bi,bj,5)          HEFFpost=SItrHEFF(i,j,bi,bj,5)
128          if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost          if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
129          SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact          SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
130       &                     +SItrFromSnow(i,j) *(1. _d 0 - growFact)       &                     +SItrFromFlood(i,j) *(1. _d 0 - growFact)
131          if ((SItrName(iTr).EQ.'age').OR.(SItrName(iTr).EQ.'one')) then  c rk: flooding can only imply an ocean-ice tracer exchange, as long
132  c Including this term in the bucket only makes sense for diagnostics purposes  c as we dont have snow tracers, so it goes through SItrBucket.
 c and should not be done for tracers that are exchanged with the ocean (we would  
 c need another bucket for ice-snow exchange, and snow tracers to start with).  
133            SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)            SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
134       &             -HEFFpost*SItrFromSnow(i,j)*(1. _d 0 - growFact)       &             -HEFFpost*SItrFromFlood(i,j)*(1. _d 0 - growFact)
135          endif  #ifdef ALLOW_SITRACER_DEBUG_DIAG
136  #ifdef ALLOW_SITRACER_DIAG          DIAGarray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr)
137          diagArray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr)       &  +SItrBucket(i,j,bi,bj,iTr)-DIAGarray(I,J,5+(iTr-1)*5)
      &  +SItrBucket(i,j,bi,bj,iTr)-diagArray(I,J,5+(iTr-1)*5)  
138  #endif  #endif
139         ENDDO         ENDDO
140        ENDDO        ENDDO
141    c TAF?      if (SItrMate(iTr).EQ.'AREA') then
142          else
143    c 1) or seaice cover expansion
144    c ============================
145    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    c ocean-ice system, the contraction/expansion terms are all
148    c simultaneous (which is sane), and the only generic effect
149    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  c 2) very ice tracer processes  c 2) very ice tracer processes
165  c ============================  c ============================
166        if (SItrName(iTr).EQ.'age') then        if (SItrName(iTr).EQ.'age') then
167  c age tracer: grow old as time passes by  c age tracer: grow old as time passes by
168         DO J=1,sNy         DO J=1,sNy
169          DO I=1,sNx          DO I=1,sNx
170            if (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0) then            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              SItracer(i,j,bi,bj,iTr)=              SItracer(i,j,bi,bj,iTr)=
174       &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm       &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
175            else            else
# Line 158  c age tracer: grow old as time passes by Line 181  c age tracer: grow old as time passes by
181  c salinity tracer: no specific process  c salinity tracer: no specific process
182        elseif (SItrName(iTr).EQ.'one') then        elseif (SItrName(iTr).EQ.'one') then
183  c "ice concentration" tracer: no specific process  c "ice concentration" tracer: no specific process
184          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          endif
202    c 3) ice-ocean tracer exchange/mapping to external variables
203    c ==========================================================
204    #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        endif
211  c 3) ice-ocean tracer exchange  #endif
 c =============================  
212        if (SItrName(iTr).EQ.'age') then        if (SItrName(iTr).EQ.'age') then
213  c age tracer: not passed to ocean  c age tracer: not passed to ocean
214        elseif (SItrName(iTr).EQ.'salinity') then        elseif (SItrName(iTr).EQ.'salinity') then
# Line 178  c "ice concentration" tracer: not passed Line 226  c "ice concentration" tracer: not passed
226        endif        endif
227        DO J=1,sNy        DO J=1,sNy
228         DO I=1,sNx         DO I=1,sNx
229  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
230          diagArray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)          DIAGarray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)
231       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
232  #endif  #endif
233  c empty bucket  c empty bucket
234          SItrBucket(i,j,bi,bj,iTr)=0. _d 0          SItrBucket(i,j,bi,bj,iTr)=0. _d 0
235         ENDDO         ENDDO
236        ENDDO        ENDDO
237    c TAF? elseif (SItrMate(iTr).EQ.'AREA') then
238  c 4) diagnostics  c 4) diagnostics
239  c ==============  c ==============
240  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
241          if (SItrMate(iTr).EQ.'HEFF') then
242        DO J=1,sNy        DO J=1,sNy
243         DO I=1,sNx         DO I=1,sNx
244          HEFFpost=SItrHEFF(i,j,bi,bj,5)          HEFFpost=SItrHEFF(i,j,bi,bj,5)
245          DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)          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          DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*HEFFpost
247  c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)  c DIAGarray(:,:,3) is the term of comparison for DIAGarray(:,:,2)
248          if (SItrName(iTr).EQ.'age') then          if (SItrName(iTr).EQ.'age') then
249            DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2)            DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2)
250          elseif (SItrName(iTr).EQ.'salinity') then          elseif (SItrName(iTr).EQ.'salinity') then
# Line 202  c diagArray(:,:,3) is the term of compar Line 252  c diagArray(:,:,3) is the term of compar
252          elseif (SItrName(iTr).EQ.'one') then          elseif (SItrName(iTr).EQ.'one') then
253            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
254          endif          endif
255  c diagArray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys  c DIAGarray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys
256  c diagArray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer)  c DIAGarray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer)
257         ENDDO         ENDDO
258        ENDDO        ENDDO
259          else
260          DO J=1,sNy
261           DO I=1,sNx
262            AREApost=SItrAREA(i,j,bi,bj,3)
263            DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
264            DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*AREApost
265    c DIAGarray(:,:,3) is the term of comparison for DIAGarray(:,:,2)
266            if (SItrName(iTr).EQ.'age') then
267              DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,1)
268            endif
269           ENDDO
270          ENDDO
271          endif
272  #endif  #endif
273        ENDDO        ENDDO
274  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
275        CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)        CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
276  #endif  #endif
277        ENDDO        ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22