/[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.9 by gforget, Thu Dec 27 23:05:47 2012 UTC
# Line 5  C $Name$ Line 5  C $Name$
5    
6  CStartOfInterface  CStartOfInterface
7        SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid )        SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid )
8  C     /=======================================================\  C     *=======================================================*
9  C     | SUBROUTINE seaice_tracer_phys                         |  C     | SUBROUTINE seaice_tracer_phys
10  C     | o Time step SItr/SItrEFF as a result of               |  C     | o Time step SItr/SItrEFF as a result of
11  C     |   seaice thermodynamics and specific tracer physics   |  C     |   seaice thermodynamics and specific tracer physics
12  C     \=======================================================/  C     *=======================================================*
13        IMPLICIT NONE        IMPLICIT NONE
14    
15  C     === Global variables ===  C     === Global variables ===
# Line 21  C     === Global variables === Line 21  C     === Global variables ===
21  #include "SEAICE.h"  #include "SEAICE.h"
22  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
23  #include "SEAICE_TRACER.h"  #include "SEAICE_TRACER.h"
24    #ifdef ALLOW_SALT_PLUME
25    # include "SALT_PLUME.h"
26    #endif
27    
28  C     === Routine arguments ===  C     === Routine arguments ===
29  C     INPUT:  C     INPUT:
# Line 37  C     === Local variables === Line 40  C     === Local variables ===
40    
41        INTEGER iTr, jTh, I, J, bi, bj, ks        INTEGER iTr, jTh, I, J, bi, bj, ks
42        _RL SItrFromOcean  (1:sNx,1:sNy)        _RL SItrFromOcean  (1:sNx,1:sNy)
43        _RL SItrFromSnow   (1:sNx,1:sNy)        _RL SItrFromFlood   (1:sNx,1:sNy)
44        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
45          _RL SItrExpand  (1:sNx,1:sNy)
46          _RL AREAprev, AREApost, expandFact
47          CHARACTER*8   diagName
48    
49  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
50        _RL DIAGarray     (1:sNx,1:sNy,Nr)        _RL DIAGarray     (1:sNx,1:sNy,Nr)
51  #endif  #endif
52    
# Line 52  cgf eventually I will need to handle the Line 58  cgf eventually I will need to handle the
58    
59        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
60        DO bi=myBxLo(myThid),myBxHi(myThid)        DO bi=myBxLo(myThid),myBxHi(myThid)
61        DO iTr=1,SItrMaxNum        DO iTr=1,SItrNumInUse
62    
63  c 0) set ice-ocean and ice-snow exchange values        c 0) set ice-ocean and ice-snow exchange values
64  c =============================================  c =============================================
65        DO J=1,sNy        DO J=1,sNy
66         DO I=1,sNx         DO I=1,sNx
67          SItrFromOcean(i,j)=0. _d 0          SItrFromOcean(i,j)=SItrFromOcean0(iTr)
68          SItrFromSnow(i,j)=0. _d 0          SItrFromFlood(i,j)=SItrFromFlood0(iTr)
69            SItrExpand(i,j)=SItrExpand0(iTr)
70         ENDDO         ENDDO
71        ENDDO        ENDDO
72        if (SItrName(iTr).EQ.'age') then  c salinity tracer:
73  c age tracer: not passed from ocean or snow        if ( (SItrName(iTr).EQ.'salinity').AND.
74         &      (SItrFromOceanFrac(iTr).GT.ZERO) ) then
75         DO J=1,sNy         DO J=1,sNy
76          DO I=1,sNx          DO I=1,sNx
77           SItrFromOcean(i,j)=0. _d 0           SItrFromOcean(i,j)=SItrFromOceanFrac(iTr)*salt(I,j,ks,bi,bj)
78          ENDDO           SItrFromFlood(i,j)=SItrFromFloodFrac(iTr)*salt(I,j,ks,bi,bj)
        ENDDO  
       elseif (SItrName(iTr).EQ.'salinity') then  
 c salinity tracer:  
        DO J=1,sNy  
         DO I=1,sNx  
          SItrFromOcean(i,j)=SIsal0  
 #ifdef SEAICE_VARIABLE_SALINITY  
          if (SIsalFRAC.GT.0.)  
      &   SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj)  
 #endif  
         ENDDO  
        ENDDO  
       elseif (SItrName(iTr).EQ.'one') then  
 c "ice concentration" tracer that should remain .EQ.1.  
        DO J=1,sNy  
         DO I=1,sNx  
          SItrFromOcean(i,j)=1. _d 0  
          SItrFromSnow(i,j)=1. _d 0  
79          ENDDO          ENDDO
80         ENDDO         ENDDO
81        endif        endif
82  c 1) seaice thermodynamics processes  c 1) seaice thermodynamics processes
83  c ==================================  c ==================================
84          if (SItrMate(iTr).EQ.'HEFF') then
85        DO J=1,sNy        DO J=1,sNy
86         DO I=1,sNx         DO I=1,sNx
87          HEFFprev=SItrHEFF(i,j,bi,bj,1)          HEFFprev=SItrHEFF(i,j,bi,bj,1)
88  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
89          diagArray(I,J,5+(iTr-1)*5) =          DIAGarray(I,J,5+(iTr-1)*5) =
90       &    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)
91  #endif  #endif
92  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 117  c apply flooding term
117          HEFFpost=SItrHEFF(i,j,bi,bj,5)          HEFFpost=SItrHEFF(i,j,bi,bj,5)
118          if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost          if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
119          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
120       &                     +SItrFromSnow(i,j) *(1. _d 0 - growFact)       &                     +SItrFromFlood(i,j) *(1. _d 0 - growFact)
121          if ((SItrName(iTr).EQ.'age').OR.(SItrName(iTr).EQ.'one')) then  c rk: flooding can only imply an ocean-ice tracer exchange, as long
122  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).  
123            SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)            SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
124       &             -HEFFpost*SItrFromSnow(i,j)*(1. _d 0 - growFact)       &             -HEFFpost*SItrFromFlood(i,j)*(1. _d 0 - growFact)
125          endif  #ifdef ALLOW_SITRACER_DEBUG_DIAG
126  #ifdef ALLOW_SITRACER_DIAG          DIAGarray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr)
127          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)  
128  #endif  #endif
129         ENDDO         ENDDO
130        ENDDO        ENDDO
131    c TAF?      if (SItrMate(iTr).EQ.'AREA') then
132          else
133    c 1) or seaice cover expansion
134    c ============================
135    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    c ocean-ice system, the contraction/expansion terms are all
138    c simultaneous (which is sane), and the only generic effect
139    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  c 2) very ice tracer processes  c 2) very ice tracer processes
155  c ============================  c ============================
156        if (SItrName(iTr).EQ.'age') then        if (SItrName(iTr).EQ.'age') then
157  c age tracer: grow old as time passes by  c age tracer: grow old as time passes by
158         DO J=1,sNy         DO J=1,sNy
159          DO I=1,sNx          DO I=1,sNx
160            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)
161         &     .EQ.'HEFF') ).OR.( (SItrAREA(i,j,bi,bj,3).GT.0. _d 0).AND.
162         &     (SItrMate(iTr).EQ.'AREA') )) then
163              SItracer(i,j,bi,bj,iTr)=              SItracer(i,j,bi,bj,iTr)=
164       &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm       &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
165            else            else
# Line 158  c age tracer: grow old as time passes by Line 171  c age tracer: grow old as time passes by
171  c salinity tracer: no specific process  c salinity tracer: no specific process
172        elseif (SItrName(iTr).EQ.'one') then        elseif (SItrName(iTr).EQ.'one') then
173  c "ice concentration" tracer: no specific process  c "ice concentration" tracer: no specific process
174          elseif (SItrName(iTr).EQ.'ridge') then
175    c simple, made up, ice surface roughness index prototype
176           DO J=1,sNy
177            DO I=1,sNx
178    c ridging increases roughness
179              SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)+
180         &    MAX(0. _d 0, SItrAREA(i,j,bi,bj,1)-SItrAREA(i,j,bi,bj,2))
181    c ice melt reduces ridges/roughness
182              HEFFprev=SItrHEFF(i,j,bi,bj,1)
183              HEFFpost=SItrHEFF(i,j,bi,bj,4)
184              tmpscal1=1. _d 0
185              if (HEFFprev.GT.HEFFpost) tmpscal1=HEFFpost/HEFFprev
186              SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*tmpscal1
187            ENDDO
188           ENDDO
189        endif        endif
190  c 3) ice-ocean tracer exchange  c 3) ice-ocean tracer exchange/mapping to external variables
191  c =============================  c ==========================================================
192        if (SItrName(iTr).EQ.'age') then  #ifdef ALLOW_DIAGNOSTICS
193  c age tracer: not passed to ocean        if (SItrMate(iTr).EQ.'HEFF') then
194        elseif (SItrName(iTr).EQ.'salinity') then          WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'Fx'
195            tmpscal1=-ONE/SEAICE_deltaTtherm*SEAICE_rhoIce
196            CALL DIAGNOSTICS_SCALE_FILL(SItrBucket(1-OLx,1-OLy,bi,bj,iTr),
197         &   tmpscal1, 1, diagName,0,1,2,bi,bj,myThid)
198          endif
199    #endif
200    
201          if ( (SItrName(iTr).EQ.'salinity').AND.
202         &     (SEAICE_salinityTracer) ) then
203  c salinity tracer: salt flux  c salinity tracer: salt flux
204  c      DO J=1,sNy          DO J=1,sNy
205  c       DO I=1,sNx           DO I=1,sNx
206  c        saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)            saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
207  c    &     *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce       &      *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
208  c note: at this point of the time step, that is the correct sign  c note: at this point of the time step, that is the correct sign
209  c        saltPlumeFlux(I,J,bi,bj) = ...  #ifdef ALLOW_SALT_PLUME
210  c       ENDDO  c should work for both constant and variable ice salinity -- to be tested
211  c      ENDDO            saltPlumeFlux(I,J,bi,bj) = MAX(ZERO,saltFlux(I,J,bi,bj))
212        elseif (SItrName(iTr).EQ.'one') then       &      *SPsalFRAC*(salt(I,j,ks,bi,bj)-SItrFromOcean(i,j))
213  c "ice concentration" tracer: not passed to ocean  #endif
214             ENDDO
215            ENDDO
216        endif        endif
217    
218        DO J=1,sNy        DO J=1,sNy
219         DO I=1,sNx         DO I=1,sNx
220  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
221          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)
222       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
223  #endif  #endif
224  c empty bucket  c empty bucket
225          SItrBucket(i,j,bi,bj,iTr)=0. _d 0          SItrBucket(i,j,bi,bj,iTr)=0. _d 0
226         ENDDO         ENDDO
227        ENDDO        ENDDO
228    
229    c TAF? elseif (SItrMate(iTr).EQ.'AREA') then
230    
231  c 4) diagnostics  c 4) diagnostics
232  c ==============  c ==============
233  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
234          if (SItrMate(iTr).EQ.'HEFF') then
235        DO J=1,sNy        DO J=1,sNy
236         DO I=1,sNx         DO I=1,sNx
237          HEFFpost=SItrHEFF(i,j,bi,bj,5)          HEFFpost=SItrHEFF(i,j,bi,bj,5)
238          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)
239          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
240  c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)  c DIAGarray(:,:,3) is the term of comparison for DIAGarray(:,:,2)
241          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  
242            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
243          elseif (SItrName(iTr).EQ.'one') then          elseif (SItrName(iTr).EQ.'one') then
244            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
245          endif          endif
246  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
247  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)
248           ENDDO
249          ENDDO
250          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         ENDDO         ENDDO
257        ENDDO        ENDDO
258          endif
259  #endif  #endif
260        ENDDO        ENDDO
261  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
262        CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)  c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
263  #endif  #endif
264        ENDDO        ENDDO
265        ENDDO        ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22