/[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.5 by gforget, Fri Feb 3 13:34:32 2012 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 54  cgf eventually I will need to handle the Line 57  cgf eventually I will need to handle the
57        DO bi=myBxLo(myThid),myBxHi(myThid)        DO bi=myBxLo(myThid),myBxHi(myThid)
58        DO iTr=1,SItrMaxNum        DO iTr=1,SItrMaxNum
59    
60  c 0) set ice-ocean and ice-snow exchange values        c 0) set ice-ocean and ice-snow exchange values
61  c =============================================  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
74          DO I=1,sNx          DO I=1,sNx
75           SItrFromOcean(i,j)=SIsal0           SItrFromOcean(i,j)=SIsal0
76  #ifdef SEAICE_VARIABLE_SALINITY  #ifdef SEAICE_VARIABLE_SALINITY
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.'salinity') then
           DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2)  
         elseif (SItrName(iTr).EQ.'salinity') then  
249            DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce            DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce
250          elseif (SItrName(iTr).EQ.'one') then          elseif (SItrName(iTr).EQ.'one') then
251            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
252          endif          endif
253  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
254  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)
255         ENDDO         ENDDO
256        ENDDO        ENDDO
257          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  #endif  #endif
267        ENDDO        ENDDO
268  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
269        CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)  c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
270  #endif  #endif
271        ENDDO        ENDDO
272        ENDDO        ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22