/[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.11 by torge, Sat Aug 23 20:22:16 2014 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 ===
16  #include "SIZE.h"  #include "SIZE.h"
17  #include "EEPARAMS.h"  #include "EEPARAMS.h"
18    #include "PARAMS.h"
19  #include "FFIELDS.h"  #include "FFIELDS.h"
20  #include "DYNVARS.h"  #include "DYNVARS.h"
21  #include "SEAICE_SIZE.h"  #include "SEAICE_SIZE.h"
22  #include "SEAICE.h"  #include "SEAICE.h"
23  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
24  #include "SEAICE_TRACER.h"  #include "SEAICE_TRACER.h"
25    #ifdef ALLOW_SALT_PLUME
26    # include "SALT_PLUME.h"
27    #endif
28    
29  C     === Routine arguments ===  C     === Routine arguments ===
30  C     INPUT:  C     INPUT:
# Line 37  C     === Local variables === Line 41  C     === Local variables ===
41    
42        INTEGER iTr, jTh, I, J, bi, bj, ks        INTEGER iTr, jTh, I, J, bi, bj, ks
43        _RL SItrFromOcean  (1:sNx,1:sNy)        _RL SItrFromOcean  (1:sNx,1:sNy)
44        _RL SItrFromSnow   (1:sNx,1:sNy)        _RL SItrFromFlood   (1:sNx,1:sNy)
45        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
46          _RL SItrExpand  (1:sNx,1:sNy)
47          _RL AREAprev, AREApost, expandFact
48          CHARACTER*8   diagName
49    
50  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
51        _RL DIAGarray     (1:sNx,1:sNy,Nr)        _RL DIAGarray     (1:sNx,1:sNy,Nr)
52  #endif  #endif
53    
# Line 52  cgf eventually I will need to handle the Line 59  cgf eventually I will need to handle the
59    
60        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
61        DO bi=myBxLo(myThid),myBxHi(myThid)        DO bi=myBxLo(myThid),myBxHi(myThid)
62        DO iTr=1,SItrMaxNum        DO iTr=1,SItrNumInUse
63    
64  c 0) set ice-ocean and ice-snow exchange values        c 0) set ice-ocean and ice-snow exchange values
65  c =============================================  c =============================================
66        DO J=1,sNy        DO J=1,sNy
67         DO I=1,sNx         DO I=1,sNx
68          SItrFromOcean(i,j)=0. _d 0          SItrFromOcean(i,j)=SItrFromOcean0(iTr)
69          SItrFromSnow(i,j)=0. _d 0          SItrFromFlood(i,j)=SItrFromFlood0(iTr)
70            SItrExpand(i,j)=SItrExpand0(iTr)
71         ENDDO         ENDDO
72        ENDDO        ENDDO
73        if (SItrName(iTr).EQ.'age') then  c salinity tracer:
74  c age tracer: not passed from ocean or snow        if ( (SItrName(iTr).EQ.'salinity').AND.
75         &      (SItrFromOceanFrac(iTr).GT.ZERO) ) then
76         DO J=1,sNy         DO J=1,sNy
77          DO I=1,sNx          DO I=1,sNx
78           SItrFromOcean(i,j)=0. _d 0           SItrFromOcean(i,j)=SItrFromOceanFrac(iTr)*salt(I,j,ks,bi,bj)
79          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  
80          ENDDO          ENDDO
81         ENDDO         ENDDO
82        endif        endif
83  c 1) seaice thermodynamics processes  c 1) seaice thermodynamics processes
84  c ==================================  c ==================================
85          if (SItrMate(iTr).EQ.'HEFF') then
86        DO J=1,sNy        DO J=1,sNy
87         DO I=1,sNx         DO I=1,sNx
88          HEFFprev=SItrHEFF(i,j,bi,bj,1)          HEFFprev=SItrHEFF(i,j,bi,bj,1)
89  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
90          diagArray(I,J,5+(iTr-1)*5) =          DIAGarray(I,J,5+(iTr-1)*5) =
91       &    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)
92  #endif  #endif
93  c apply the sequence of thermodynamics increments to actual traceur  c apply the sequence of thermodynamics increments to actual tracer
94  c (see seaice_growth.F)  c (see seaice_growth.F)
95  c (jTh=1 tendency due to ice-ocean interaction)  c (jTh=1 tendency due to ice-ocean interaction)
96  c (jTh=2 tendency due to the atmosphere, over ice covered part)  c (jTh=2 tendency due to the atmosphere, over ice covered part)
# Line 126  c apply flooding term Line 118  c apply flooding term
118          HEFFpost=SItrHEFF(i,j,bi,bj,5)          HEFFpost=SItrHEFF(i,j,bi,bj,5)
119          if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost          if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
120          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
121       &                     +SItrFromSnow(i,j) *(1. _d 0 - growFact)       &                     +SItrFromFlood(i,j) *(1. _d 0 - growFact)
122          if ((SItrName(iTr).EQ.'age').OR.(SItrName(iTr).EQ.'one')) then  c rk: flooding can only imply an ocean-ice tracer exchange, as long
123  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).  
124            SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)            SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
125       &             -HEFFpost*SItrFromSnow(i,j)*(1. _d 0 - growFact)       &             -HEFFpost*SItrFromFlood(i,j)*(1. _d 0 - growFact)
126          endif  #ifdef ALLOW_SITRACER_DEBUG_DIAG
127  #ifdef ALLOW_SITRACER_DIAG          DIAGarray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr)
128          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)  
129  #endif  #endif
130         ENDDO         ENDDO
131        ENDDO        ENDDO
132    c TAF?      if (SItrMate(iTr).EQ.'AREA') then
133          else
134    c 1) or seaice cover expansion
135    c ============================
136    c this is much simpler than for ice volume/mass tracers, because
137    c properties of the ice surface are not be conserved across the
138    c ocean-ice system, the contraction/expansion terms are all
139    c simultaneous (which is sane), and the only generic effect
140    c is due to expansion (new cover).
141          DO J=1,sNy
142           DO I=1,sNx
143    c apply expansion
144            AREAprev=SItrAREA(i,j,bi,bj,2)
145            AREApost=SItrAREA(i,j,bi,bj,3)
146    c compute ratio in [0. 1.] range for expansion/contraction
147            expandFact=1. _d 0
148            if (AREApost.GT.AREAprev) expandFact=AREAprev/AREApost
149    c update SItr accordingly
150             SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*expandFact
151         &                      +SItrExpand(i,j)*(1. _d 0 - expandFact)
152           ENDDO
153          ENDDO
154          endif
155  c 2) very ice tracer processes  c 2) very ice tracer processes
156  c ============================  c ============================
157        if (SItrName(iTr).EQ.'age') then        if (SItrName(iTr).EQ.'age') then
158  c age tracer: grow old as time passes by  c age tracer: grow old as time passes by
159         DO J=1,sNy         DO J=1,sNy
160          DO I=1,sNx          DO I=1,sNx
161            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)
162         &     .EQ.'HEFF') ).OR.( (SItrAREA(i,j,bi,bj,3).GT.0. _d 0).AND.
163         &     (SItrMate(iTr).EQ.'AREA') )) then
164              SItracer(i,j,bi,bj,iTr)=              SItracer(i,j,bi,bj,iTr)=
165       &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm       &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
166            else            else
# Line 158  c age tracer: grow old as time passes by Line 172  c age tracer: grow old as time passes by
172  c salinity tracer: no specific process  c salinity tracer: no specific process
173        elseif (SItrName(iTr).EQ.'one') then        elseif (SItrName(iTr).EQ.'one') then
174  c "ice concentration" tracer: no specific process  c "ice concentration" tracer: no specific process
175          elseif (SItrName(iTr).EQ.'ridge') then
176    c simple, made up, ice surface roughness index prototype
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        endif
191  c 3) ice-ocean tracer exchange  c 3) ice-ocean tracer exchange/mapping to external variables
192  c =============================  c ==========================================================
193        if (SItrName(iTr).EQ.'age') then  #ifdef ALLOW_DIAGNOSTICS
194  c age tracer: not passed to ocean        IF ( useDiagnostics .AND. SItrMate(iTr).EQ.'HEFF') THEN
195        elseif (SItrName(iTr).EQ.'salinity') then          WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'Fx'
196            tmpscal1=-ONE/SEAICE_deltaTtherm*SEAICE_rhoIce
197            CALL DIAGNOSTICS_SCALE_FILL(SItrBucket(1-OLx,1-OLy,bi,bj,iTr),
198         &   tmpscal1, 1, diagName,0,1,2,bi,bj,myThid)
199          ENDIF
200    #endif
201    
202          if ( (SItrName(iTr).EQ.'salinity').AND.
203         &     (SEAICE_salinityTracer) ) then
204  c salinity tracer: salt flux  c salinity tracer: salt flux
205  c      DO J=1,sNy          DO J=1,sNy
206  c       DO I=1,sNx           DO I=1,sNx
207  c        saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)            saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
208  c    &     *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce       &      *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
209  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
210  c        saltPlumeFlux(I,J,bi,bj) = ...  #ifdef ALLOW_SALT_PLUME
211  c       ENDDO  c should work for both constant and variable ice salinity -- to be tested
212  c      ENDDO            saltPlumeFlux(I,J,bi,bj) = MAX(ZERO,saltFlux(I,J,bi,bj))
213        elseif (SItrName(iTr).EQ.'one') then       &      *SPsalFRAC*(salt(I,j,ks,bi,bj)-SItrFromOcean(i,j))
214  c "ice concentration" tracer: not passed to ocean  #endif
215             ENDDO
216            ENDDO
217        endif        endif
218    
219        DO J=1,sNy        DO J=1,sNy
220         DO I=1,sNx         DO I=1,sNx
221  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
222          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)
223       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
224  #endif  #endif
225  c empty bucket  c empty bucket
226          SItrBucket(i,j,bi,bj,iTr)=0. _d 0  c  but not for 'grease' (see seaice_growth.F)
227            if (SItrName(iTr).NE.'grease')
228         &      SItrBucket(i,j,bi,bj,iTr)=0. _d 0
229         ENDDO         ENDDO
230        ENDDO        ENDDO
231    
232    c TAF? elseif (SItrMate(iTr).EQ.'AREA') then
233    
234  c 4) diagnostics  c 4) diagnostics
235  c ==============  c ==============
236  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
237          if (SItrMate(iTr).EQ.'HEFF') then
238        DO J=1,sNy        DO J=1,sNy
239         DO I=1,sNx         DO I=1,sNx
240          HEFFpost=SItrHEFF(i,j,bi,bj,5)          HEFFpost=SItrHEFF(i,j,bi,bj,5)
241          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)
242          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
243  c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)  c DIAGarray(:,:,3) is the term of comparison for DIAGarray(:,:,2)
244          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  
245            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
246          elseif (SItrName(iTr).EQ.'one') then          elseif (SItrName(iTr).EQ.'one') then
247            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
248          endif          endif
249  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
250  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)
251           ENDDO
252          ENDDO
253          else
254          DO J=1,sNy
255           DO I=1,sNx
256            AREApost=SItrAREA(i,j,bi,bj,3)
257            DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
258            DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*AREApost
259         ENDDO         ENDDO
260        ENDDO        ENDDO
261          endif
262  #endif  #endif
263        ENDDO        ENDDO
264  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
265        CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)  c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
266  #endif  #endif
267        ENDDO        ENDDO
268        ENDDO        ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22