/[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.2 by gforget, Thu Jun 9 19:37:01 2011 UTC revision 1.8 by gforget, Mon Mar 5 14:52:16 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 41  C     === Local variables === Line 44  C     === Local variables ===
44        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1        _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
45        _RL SItrExpand  (1:sNx,1:sNy)        _RL SItrExpand  (1:sNx,1:sNy)
46        _RL AREAprev, AREApost, expandFact        _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 54  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          SItrFromFlood(i,j)=0. _d 0          SItrFromFlood(i,j)=SItrFromFlood0(iTr)
69          SItrExpand(i,j)=0. _d 0          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: no age in ocean, or effect from ice cover changes        if ( (SItrName(iTr).EQ.'salinity').AND.
74        elseif (SItrName(iTr).EQ.'salinity') then       &      (SItrFromOceanFrac(iTr).GT.ZERO) ) then
 c salinity tracer:  
75         DO J=1,sNy         DO J=1,sNy
76          DO I=1,sNx          DO I=1,sNx
77           SItrFromOcean(i,j)=SIsal0           SItrFromOcean(i,j)=SItrFromOceanFrac(iTr)*salt(I,j,ks,bi,bj)
78  #ifdef SEAICE_VARIABLE_SALINITY           SItrFromFlood(i,j)=SItrFromFloodFrac(iTr)*salt(I,j,ks,bi,bj)
          if (SIsalFRAC.GT.0.)  
      &   SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj)  
 #endif  
 c as of now, flooding implies no salt extraction from ocean  
         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  
          SItrFromFlood(i,j)=1. _d 0  
79          ENDDO          ENDDO
80         ENDDO         ENDDO
81        endif        endif
# Line 94  c ================================== Line 85  c ==================================
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 131  c rk: flooding can only imply an ocean-i Line 122  c rk: flooding can only imply an ocean-i
122  c as we dont have snow tracers, so it goes through SItrBucket.  c as we dont have snow tracers, so it goes through SItrBucket.
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*SItrFromFlood(i,j)*(1. _d 0 - growFact)       &             -HEFFpost*SItrFromFlood(i,j)*(1. _d 0 - growFact)
125  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
126          diagArray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr)          DIAGarray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr)
127       &  +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
# Line 141  c TAF?      if (SItrMate(iTr).EQ.'AREA') Line 132  c TAF?      if (SItrMate(iTr).EQ.'AREA')
132        else        else
133  c 1) or seaice cover expansion  c 1) or seaice cover expansion
134  c ============================  c ============================
135  c this is much simpler than for ice volume/mass tracers, because  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  c properties of the ice surface are not be conserved across the
137  c ocean-ice system, the contraction/expansion terms are all  c ocean-ice system, the contraction/expansion terms are all
138  c simultaneous (which is sane), and the only generic effect  c simultaneous (which is sane), and the only generic effect
139  c is due to expansion (new cover).  c is due to expansion (new cover).
140        DO J=1,sNy        DO J=1,sNy
141         DO I=1,sNx         DO I=1,sNx
# Line 200  c ice melt reduces ridges/roughness Line 191  c ice melt reduces ridges/roughness
191        endif        endif
192  c 3) ice-ocean tracer exchange/mapping to external variables  c 3) ice-ocean tracer exchange/mapping to external variables
193  c ==========================================================  c ==========================================================
194        if (SItrName(iTr).EQ.'age') then  #ifdef ALLOW_DIAGNOSTICS
195  c age tracer: not passed to ocean        if (SItrMate(iTr).EQ.'HEFF') then
196        elseif (SItrName(iTr).EQ.'salinity') then          WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'Fx'
197            tmpscal1=-ONE/SEAICE_deltaTtherm*SEAICE_rhoIce
198            CALL DIAGNOSTICS_SCALE_FILL(SItrBucket(1-OLx,1-OLy,bi,bj,iTr),
199         &   tmpscal1, 1, diagName,0,1,2,bi,bj,myThid)
200          endif
201    #endif
202    
203          if ( (SItrName(iTr).EQ.'salinity').AND.
204         &     (SEAICE_salinityTracer) ) then
205  c salinity tracer: salt flux  c salinity tracer: salt flux
206  c      DO J=1,sNy          DO J=1,sNy
207  c       DO I=1,sNx           DO I=1,sNx
208  c        saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)            saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
209  c    &     *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce       &      *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
210  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
211  c        saltPlumeFlux(I,J,bi,bj) = ...  #ifdef ALLOW_SALT_PLUME
212  c       ENDDO  c should work for both constant and variable ice salinity -- to be tested
213  c      ENDDO            saltPlumeFlux(I,J,bi,bj) = MAX(ZERO,saltFlux(I,J,bi,bj))
214        elseif (SItrName(iTr).EQ.'one') then       &      *SPsalFRAC*(salt(I,j,ks,bi,bj)-SItrFromOcean(i,j))
215  c "ice concentration" tracer: not passed to ocean  #endif
216             ENDDO
217            ENDDO
218        endif        endif
219    
220        DO J=1,sNy        DO J=1,sNy
221         DO I=1,sNx         DO I=1,sNx
222  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
223          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)
224       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce       &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
225  #endif  #endif
226  c empty bucket  c empty bucket
227          SItrBucket(i,j,bi,bj,iTr)=0. _d 0          SItrBucket(i,j,bi,bj,iTr)=0. _d 0
228         ENDDO         ENDDO
229        ENDDO        ENDDO
230    
231  c TAF? elseif (SItrMate(iTr).EQ.'AREA') then  c TAF? elseif (SItrMate(iTr).EQ.'AREA') then
232    
233  c 4) diagnostics  c 4) diagnostics
234  c ==============  c ==============
235  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
236        if (SItrMate(iTr).EQ.'HEFF') then        if (SItrMate(iTr).EQ.'HEFF') then
237        DO J=1,sNy        DO J=1,sNy
238         DO I=1,sNx         DO I=1,sNx
239          HEFFpost=SItrHEFF(i,j,bi,bj,5)          HEFFpost=SItrHEFF(i,j,bi,bj,5)
240          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)
241          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
242  c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)  c DIAGarray(:,:,3) is the term of comparison for DIAGarray(:,:,2)
243          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  
244            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
245          elseif (SItrName(iTr).EQ.'one') then          elseif (SItrName(iTr).EQ.'one') then
246            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost            DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
247          endif          endif
248  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
249  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)
250         ENDDO         ENDDO
251        ENDDO        ENDDO
252        else        else
# Line 253  c diagArray(:,:,5) is the tracer flux fr Line 255  c diagArray(:,:,5) is the tracer flux fr
255          AREApost=SItrAREA(i,j,bi,bj,3)          AREApost=SItrAREA(i,j,bi,bj,3)
256          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)
257          DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*AREApost          DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*AREApost
 c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)  
         if (SItrName(iTr).EQ.'age') then  
           DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,1)  
         endif  
258         ENDDO         ENDDO
259        ENDDO        ENDDO
260        endif        endif
261  #endif  #endif
262        ENDDO        ENDDO
263  #ifdef ALLOW_SITRACER_DIAG  #ifdef ALLOW_SITRACER_DEBUG_DIAG
264        CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)  c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
265  #endif  #endif
266        ENDDO        ENDDO
267        ENDDO        ENDDO
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22