pkg/seaice/seaice_tracer_phys.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_tracer_growth.F,v 1.1 2011/06/02 16:21:08 gforget Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid )
C     /=======================================================\
C     | SUBROUTINE seaice_tracer_phys                         |
C     | o Time step SItr/SItrEFF as a result of               |
C     |   seaice thermodynamics and specific tracer physics   |
C     \=======================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "FFIELDS.h"
#include "DYNVARS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE_TRACER.h"

C     === Routine arguments ===
C     INPUT:
C     myTime  :: Simulation time
C     myIter  :: Simulation timestep number
C     myThid  :: Thread no. that called this routine.
C     OUTPUT:
      _RL myTime
      INTEGER myIter, myThid
CEndOfInterface

C     === Local variables ===
#ifdef ALLOW_SITRACER

      INTEGER iTr, jTh, I, J, bi, bj, ks
      _RL SItrFromOcean  (1:sNx,1:sNy)
      _RL SItrFromSnow   (1:sNx,1:sNy)
      _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1

#ifdef ALLOW_SITRACER_DIAG
      _RL DIAGarray     (1:sNx,1:sNy,Nr)
#endif

cgf for now I do not fully account for ocean-ice fluxes of tracer
cgf -> I just prescribe it consistent with age tracer
cgf eventually I will need to handle them as function params

      ks=1

      DO bj=myByLo(myThid),myByHi(myThid)
      DO bi=myBxLo(myThid),myBxHi(myThid)
      DO iTr=1,SItrMaxNum

c 0) set ice-ocean and ice-snow exchange values       
c =============================================
      DO J=1,sNy
       DO I=1,sNx
        SItrFromOcean(i,j)=0. _d 0
        SItrFromSnow(i,j)=0. _d 0
       ENDDO
      ENDDO
      if (SItrName(iTr).EQ.'age') then
c age tracer: not passed from ocean or snow
       DO J=1,sNy
        DO I=1,sNx
         SItrFromOcean(i,j)=0. _d 0
        ENDDO
       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
        ENDDO
       ENDDO
      endif
c 1) seaice thermodynamics processes
c ==================================
      DO J=1,sNy
       DO I=1,sNx
        HEFFprev=SItrHEFF(i,j,bi,bj,1)
#ifdef ALLOW_SITRACER_DIAG
        diagArray(I,J,5+(iTr-1)*5) = 
     &    HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr)
#endif
c apply the sequence of thermodynamics increments to actual traceur
c (see seaice_growth.F)
c (jTh=1 tendency due to ice-ocean interaction)
c (jTh=2 tendency due to the atmosphere, over ice covered part)
c (jTh=3 tendency due to the atmosphere, over open water part)
c (jTh=4 tendency due to flooding)
        DO jTh=1,3
         HEFFprev=SItrHEFF(i,j,bi,bj,jTh)
         HEFFpost=SItrHEFF(i,j,bi,bj,jTh+1)
c compute ratio in [0. 1.] range for either growth or melt
         growFact=1. _d 0
         meltPart=0. _d 0
         if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
         if (HEFFpost.LT.HEFFprev) meltPart=HEFFprev-HEFFpost
c update SItr accordingly
         SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
     &                      +SItrFromOcean(i,j)*(1. _d 0 - growFact)
         SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &             -HEFFpost*SItrFromOcean(i,j)*(1. _d 0 - growFact)
         SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &                            +meltPart*SItracer(i,j,bi,bj,iTr)
        ENDDO
c apply flooding term
        growFact=1. _d 0
        HEFFprev=SItrHEFF(i,j,bi,bj,4)
        HEFFpost=SItrHEFF(i,j,bi,bj,5)
        if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
        SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
     &                     +SItrFromSnow(i,j) *(1. _d 0 - growFact)
        if ((SItrName(iTr).EQ.'age').OR.(SItrName(iTr).EQ.'one')) then
c Including this term in the bucket only makes sense for diagnostics purposes
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).
          SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &             -HEFFpost*SItrFromSnow(i,j)*(1. _d 0 - growFact)
        endif
#ifdef ALLOW_SITRACER_DIAG
        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)
#endif
       ENDDO
      ENDDO
c 2) very ice tracer processes
c ============================
      if (SItrName(iTr).EQ.'age') then
c age tracer: grow old as time passes by
       DO J=1,sNy
        DO I=1,sNx
          if (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0) then
            SItracer(i,j,bi,bj,iTr)=
     &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
          else
            SItracer(i,j,bi,bj,iTr)=0. _d 0
          endif
        ENDDO
       ENDDO
      elseif (SItrName(iTr).EQ.'salinity') then
c salinity tracer: no specific process
      elseif (SItrName(iTr).EQ.'one') then
c "ice concentration" tracer: no specific process
      endif
c 3) ice-ocean tracer exchange
c =============================
      if (SItrName(iTr).EQ.'age') then
c age tracer: not passed to ocean
      elseif (SItrName(iTr).EQ.'salinity') then
c salinity tracer: salt flux
c      DO J=1,sNy
c       DO I=1,sNx
c        saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
c    &     *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
c note: at this point of the time step, that is the correct sign
c        saltPlumeFlux(I,J,bi,bj) = ...
c       ENDDO
c      ENDDO
      elseif (SItrName(iTr).EQ.'one') then
c "ice concentration" tracer: not passed to ocean
      endif
      DO J=1,sNy
       DO I=1,sNx
#ifdef ALLOW_SITRACER_DIAG
        diagArray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)
     &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
#endif
c empty bucket
        SItrBucket(i,j,bi,bj,iTr)=0. _d 0
       ENDDO
      ENDDO
c 4) diagnostics
c ==============
#ifdef ALLOW_SITRACER_DIAG
      DO J=1,sNy
       DO I=1,sNx
        HEFFpost=SItrHEFF(i,j,bi,bj,5)
        DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
        DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*HEFFpost
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,2)
        elseif (SItrName(iTr).EQ.'salinity') then
          DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce
        elseif (SItrName(iTr).EQ.'one') then
          DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
        endif
c diagArray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys
c diagArray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer)
       ENDDO
      ENDDO
#endif
      ENDDO
#ifdef ALLOW_SITRACER_DIAG
      CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
#endif
      ENDDO
      ENDDO

#endif /* ALLOW_SITRACER */

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_tracer_growth.F,v 1.1 2011/06/02 16:21:08 gforget Exp $
2	C $Name: $
3
4	#include "SEAICE_OPTIONS.h"
5
6	CStartOfInterface
7	SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid )
8	C /=======================================================\
9	C \| SUBROUTINE seaice_tracer_phys \|
10	C \| o Time step SItr/SItrEFF as a result of \|
11	C \| seaice thermodynamics and specific tracer physics \|
12	C \=======================================================/
13	IMPLICIT NONE
14
15	C === Global variables ===
16	#include "SIZE.h"
17	#include "EEPARAMS.h"
18	#include "FFIELDS.h"
19	#include "DYNVARS.h"
20	#include "SEAICE_SIZE.h"
21	#include "SEAICE.h"
22	#include "SEAICE_PARAMS.h"
23	#include "SEAICE_TRACER.h"
24
25	C === Routine arguments ===
26	C INPUT:
27	C myTime :: Simulation time
28	C myIter :: Simulation timestep number
29	C myThid :: Thread no. that called this routine.
30	C OUTPUT:
31	_RL myTime
32	INTEGER myIter, myThid
33	CEndOfInterface
34
35	C === Local variables ===
36	#ifdef ALLOW_SITRACER
37
38	INTEGER iTr, jTh, I, J, bi, bj, ks
39	_RL SItrFromOcean (1:sNx,1:sNy)
40	_RL SItrFromSnow (1:sNx,1:sNy)
41	_RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
42
43	#ifdef ALLOW_SITRACER_DIAG
44	_RL DIAGarray (1:sNx,1:sNy,Nr)
45	#endif
46
47	cgf for now I do not fully account for ocean-ice fluxes of tracer
48	cgf -> I just prescribe it consistent with age tracer
49	cgf eventually I will need to handle them as function params
50
51	ks=1
52
53	DO bj=myByLo(myThid),myByHi(myThid)
54	DO bi=myBxLo(myThid),myBxHi(myThid)
55	DO iTr=1,SItrMaxNum
56
57	c 0) set ice-ocean and ice-snow exchange values
58	c =============================================
59	DO J=1,sNy
60	DO I=1,sNx
61	SItrFromOcean(i,j)=0. _d 0
62	SItrFromSnow(i,j)=0. _d 0
63	ENDDO
64	ENDDO
65	if (SItrName(iTr).EQ.'age') then
66	c age tracer: not passed from ocean or snow
67	DO J=1,sNy
68	DO I=1,sNx
69	SItrFromOcean(i,j)=0. _d 0
70	ENDDO
71	ENDDO
72	elseif (SItrName(iTr).EQ.'salinity') then
73	c salinity tracer:
74	DO J=1,sNy
75	DO I=1,sNx
76	SItrFromOcean(i,j)=SIsal0
77	#ifdef SEAICE_VARIABLE_SALINITY
78	if (SIsalFRAC.GT.0.)
79	& SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj)
80	#endif
81	ENDDO
82	ENDDO
83	elseif (SItrName(iTr).EQ.'one') then
84	c "ice concentration" tracer that should remain .EQ.1.
85	DO J=1,sNy
86	DO I=1,sNx
87	SItrFromOcean(i,j)=1. _d 0
88	SItrFromSnow(i,j)=1. _d 0
89	ENDDO
90	ENDDO
91	endif
92	c 1) seaice thermodynamics processes
93	c ==================================
94	DO J=1,sNy
95	DO I=1,sNx
96	HEFFprev=SItrHEFF(i,j,bi,bj,1)
97	#ifdef ALLOW_SITRACER_DIAG
98	diagArray(I,J,5+(iTr-1)*5) =
99	& HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr)
100	#endif
101	c apply the sequence of thermodynamics increments to actual traceur
102	c (see seaice_growth.F)
103	c (jTh=1 tendency due to ice-ocean interaction)
104	c (jTh=2 tendency due to the atmosphere, over ice covered part)
105	c (jTh=3 tendency due to the atmosphere, over open water part)
106	c (jTh=4 tendency due to flooding)
107	DO jTh=1,3
108	HEFFprev=SItrHEFF(i,j,bi,bj,jTh)
109	HEFFpost=SItrHEFF(i,j,bi,bj,jTh+1)
110	c compute ratio in [0. 1.] range for either growth or melt
111	growFact=1. _d 0
112	meltPart=0. _d 0
113	if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
114	if (HEFFpost.LT.HEFFprev) meltPart=HEFFprev-HEFFpost
115	c update SItr accordingly
116	SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
117	& +SItrFromOcean(i,j)*(1. _d 0 - growFact)
118	SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
119	& -HEFFpostSItrFromOcean(i,j)(1. _d 0 - growFact)
120	SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
121	& +meltPart*SItracer(i,j,bi,bj,iTr)
122	ENDDO
123	c apply flooding term
124	growFact=1. _d 0
125	HEFFprev=SItrHEFF(i,j,bi,bj,4)
126	HEFFpost=SItrHEFF(i,j,bi,bj,5)
127	if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
128	SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
129	& +SItrFromSnow(i,j) *(1. _d 0 - growFact)
130	if ((SItrName(iTr).EQ.'age').OR.(SItrName(iTr).EQ.'one')) then
131	c Including this term in the bucket only makes sense for diagnostics purposes
132	c and should not be done for tracers that are exchanged with the ocean (we would
133	c need another bucket for ice-snow exchange, and snow tracers to start with).
134	SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
135	& -HEFFpostSItrFromSnow(i,j)(1. _d 0 - growFact)
136	endif
137	#ifdef ALLOW_SITRACER_DIAG
138	diagArray(I,J,5+(iTr-1)5) = HEFFpostSItracer(i,j,bi,bj,iTr)
139	& +SItrBucket(i,j,bi,bj,iTr)-diagArray(I,J,5+(iTr-1)*5)
140	#endif
141	ENDDO
142	ENDDO
143	c 2) very ice tracer processes
144	c ============================
145	if (SItrName(iTr).EQ.'age') then
146	c age tracer: grow old as time passes by
147	DO J=1,sNy
148	DO I=1,sNx
149	if (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0) then
150	SItracer(i,j,bi,bj,iTr)=
151	& SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
152	else
153	SItracer(i,j,bi,bj,iTr)=0. _d 0
154	endif
155	ENDDO
156	ENDDO
157	elseif (SItrName(iTr).EQ.'salinity') then
158	c salinity tracer: no specific process
159	elseif (SItrName(iTr).EQ.'one') then
160	c "ice concentration" tracer: no specific process
161	endif
162	c 3) ice-ocean tracer exchange
163	c =============================
164	if (SItrName(iTr).EQ.'age') then
165	c age tracer: not passed to ocean
166	elseif (SItrName(iTr).EQ.'salinity') then
167	c salinity tracer: salt flux
168	c DO J=1,sNy
169	c DO I=1,sNx
170	c saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
171	c & HEFFM(I,J,bi,bj)/SEAICE_deltaTthermSEAICE_rhoIce
172	c note: at this point of the time step, that is the correct sign
173	c saltPlumeFlux(I,J,bi,bj) = ...
174	c ENDDO
175	c ENDDO
176	elseif (SItrName(iTr).EQ.'one') then
177	c "ice concentration" tracer: not passed to ocean
178	endif
179	DO J=1,sNy
180	DO I=1,sNx
181	#ifdef ALLOW_SITRACER_DIAG
182	diagArray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)
183	& HEFFM(I,J,bi,bj)/SEAICE_deltaTthermSEAICE_rhoIce
184	#endif
185	c empty bucket
186	SItrBucket(i,j,bi,bj,iTr)=0. _d 0
187	ENDDO
188	ENDDO
189	c 4) diagnostics
190	c ==============
191	#ifdef ALLOW_SITRACER_DIAG
192	DO J=1,sNy
193	DO I=1,sNx
194	HEFFpost=SItrHEFF(i,j,bi,bj,5)
195	DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
196	DIAGarray(I,J,2+(iTr-1)5) = SItracer(i,j,bi,bj,iTr)HEFFpost
197	c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)
198	if (SItrName(iTr).EQ.'age') then
199	DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2)
200	elseif (SItrName(iTr).EQ.'salinity') then
201	DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce
202	elseif (SItrName(iTr).EQ.'one') then
203	DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
204	endif
205	c diagArray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys
206	c diagArray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer)
207	ENDDO
208	ENDDO
209	#endif
210	ENDDO
211	#ifdef ALLOW_SITRACER_DIAG
212	CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1 ',0,Nr,3,bi,bj,myThid)
213	#endif
214	ENDDO
215	ENDDO
216
217	#endif /* ALLOW_SITRACER */
218
219	RETURN
220	END