/[MITgcm]/MITgcm/pkg/seaice/seaice_growth.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_growth.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.68 by dimitri, Fri Sep 3 23:33:49 2010 UTC revision 1.69 by dimitri, Tue Sep 14 15:14:07 2010 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "SEAICE_OPTIONS.h"  #include "SEAICE_OPTIONS.h"
5    
6  CBOP  CBOP
7  C     !ROUTINE: SEAICE_GROWTH  C     !ROUTINE: SEAICE_GROWTH
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )        SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
10  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
11  C     *==========================================================*  C     *==========================================================*
12  C     | SUBROUTINE seaice_growth  C     | SUBROUTINE seaice_growth
13  C     | o Updata ice thickness and snow depth  C     | o Updata ice thickness and snow depth
14  C     *==========================================================*  C     *==========================================================*
15  C     \ev  C     \ev
16    
17  C     !USES:  C     !USES:
18        IMPLICIT NONE        IMPLICIT NONE
19  C     === Global variables ===  C     === Global variables ===
20  #include "SIZE.h"  #include "SIZE.h"
21  #include "EEPARAMS.h"  #include "EEPARAMS.h"
22  #include "PARAMS.h"  #include "PARAMS.h"
23  #include "DYNVARS.h"  #include "DYNVARS.h"
24  #include "GRID.h"  #include "GRID.h"
25  #include "FFIELDS.h"  #include "FFIELDS.h"
26  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
27  #include "SEAICE.h"  #include "SEAICE.h"
28  #ifdef ALLOW_EXF  #ifdef ALLOW_EXF
29  # include "EXF_OPTIONS.h"  # include "EXF_OPTIONS.h"
30  # include "EXF_FIELDS.h"  # include "EXF_FIELDS.h"
31  # include "EXF_PARAM.h"  # include "EXF_PARAM.h"
32  #endif  #endif
33  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
34  # include "SALT_PLUME.h"  # include "SALT_PLUME.h"
35  #endif  #endif
36  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
37  # include "tamc.h"  # include "tamc.h"
38  #endif  #endif
39    
40  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
41  C     === Routine arguments ===  C     === Routine arguments ===
42  C     myTime :: Simulation time  C     myTime :: Simulation time
43  C     myIter :: Simulation timestep number  C     myIter :: Simulation timestep number
44  C     myThid :: Thread no. that called this routine.  C     myThid :: Thread no. that called this routine.
45        _RL myTime        _RL myTime
46        INTEGER myIter, myThid        INTEGER myIter, myThid
47  CEOP  CEOP
48    
49  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
50  C     === Local variables ===  C     === Local variables ===
51  C     i,j,bi,bj :: Loop counters  C     i,j,bi,bj :: Loop counters
52        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
53  C     number of surface interface layer  C     number of surface interface layer
54        INTEGER kSurface        INTEGER kSurface
55  C     constants  C     constants
56        _RL TBC, SDF, ICE2SNOW        _RL TBC, SDF, ICE2SNOW
57        _RL QI, recip_QI, QS        _RL QI, recip_QI, QS
58  C     auxillary variables  C     auxillary variables
59        _RL snowEnergy        _RL snowEnergy
60        _RL growthHEFF, growthNeg        _RL growthHEFF, growthNeg
61  #ifdef ALLOW_SEAICE_FLOODING  #ifdef ALLOW_SEAICE_FLOODING
62        _RL hDraft        _RL hDraft
63  #endif /* ALLOW_SEAICE_FLOODING */  #endif /* ALLOW_SEAICE_FLOODING */
64        _RL availHeat     (1:sNx,1:sNy)        _RL availHeat     (1:sNx,1:sNy)
65        _RL hEffOld       (1:sNx,1:sNy)        _RL hEffOld       (1:sNx,1:sNy)
66        _RL GAREA         (1:sNx,1:sNy)        _RL GAREA         (1:sNx,1:sNy)
67        _RL GHEFF         (1:sNx,1:sNy)        _RL GHEFF         (1:sNx,1:sNy)
68  C     RESID_HEAT is residual heat above freezing in equivalent m of ice  C     RESID_HEAT is residual heat above freezing in equivalent m of ice
69        _RL RESID_HEAT    (1:sNx,1:sNy)        _RL RESID_HEAT    (1:sNx,1:sNy)
70  #ifdef SEAICE_SALINITY  #ifdef SEAICE_SALINITY
71        _RL saltFluxAdjust(1:sNx,1:sNy)        _RL saltFluxAdjust(1:sNx,1:sNy)
72  #endif  #endif
73    
74  C     FICE  - thermodynamic ice growth rate over sea ice in W/m^2  C     FICE  - thermodynamic ice growth rate over sea ice in W/m^2
75  C             >0 causes ice growth, <0 causes snow and sea ice melt  C             >0 causes ice growth, <0 causes snow and sea ice melt
76  C     FHEFF - effective thermodynamic ice growth rate over sea ice in W/m^2  C     FHEFF - effective thermodynamic ice growth rate over sea ice in W/m^2
77  C             >0 causes ice growth, <0 causes snow and sea ice melt  C             >0 causes ice growth, <0 causes snow and sea ice melt
78  C     QNETO  - thermodynamic ice growth rate over open water in W/m^2  C     QNETO  - thermodynamic ice growth rate over open water in W/m^2
79  C             ( = surface heat flux )  C             ( = surface heat flux )
80  C             >0 causes ice growth, <0 causes snow and sea ice melt  C             >0 causes ice growth, <0 causes snow and sea ice melt
81  C     QNETI  - net surface heat flux under ice in W/m^2  C     QNETI  - net surface heat flux under ice in W/m^2
82  C     QSWO   - short wave heat flux over ocean in W/m^2  C     QSWO   - short wave heat flux over ocean in W/m^2
83  C     QSWI   - short wave heat flux under ice in W/m^2  C     QSWI   - short wave heat flux under ice in W/m^2
84        _RL FHEFF         (1:sNx,1:sNy)        _RL FHEFF         (1:sNx,1:sNy)
85        _RL FICE          (1:sNx,1:sNy)        _RL FICE          (1:sNx,1:sNy)
86        _RL QNETO         (1:sNx,1:sNy)        _RL QNETO         (1:sNx,1:sNy)
87        _RL QNETI         (1:sNx,1:sNy)        _RL QNETI         (1:sNx,1:sNy)
88        _RL QSWO          (1:sNx,1:sNy)        _RL QSWO          (1:sNx,1:sNy)
89        _RL QSWI          (1:sNx,1:sNy)        _RL QSWI          (1:sNx,1:sNy)
90  C  C
91        _RL HCORR         (1:sNx,1:sNy)        _RL HCORR         (1:sNx,1:sNy)
92  C     actual ice thickness with upper and lower limit  C     actual ice thickness with upper and lower limit
93        _RL HICE          (1:sNx,1:sNy)        _RL HICE          (1:sNx,1:sNy)
94  C     actual snow thickness  C     actual snow thickness
95        _RL hSnwLoc       (1:sNx,1:sNy)        _RL hSnwLoc       (1:sNx,1:sNy)
96  C     wind speed  C     wind speed
97        _RL UG            (1:sNx,1:sNy)        _RL UG            (1:sNx,1:sNy)
98        _RL SPEED_SQ        _RL SPEED_SQ
99  C     local copy of AREA  C     local copy of AREA
100        _RL areaLoc        _RL areaLoc
101    
102  #ifdef SEAICE_MULTICATEGORY  #ifdef SEAICE_MULTICATEGORY
103        INTEGER it        INTEGER it
104        INTEGER ilockey        INTEGER ilockey
105        _RL RK        _RL RK
106        _RL HICEP         (1:sNx,1:sNy)        _RL HICEP         (1:sNx,1:sNy)
107        _RL FICEP         (1:sNx,1:sNy)        _RL FICEP         (1:sNx,1:sNy)
108        _RL QSWIP         (1:sNx,1:sNy)        _RL QSWIP         (1:sNx,1:sNy)
109  #endif  #endif
110    
111  #ifdef SEAICE_AGE  #ifdef SEAICE_AGE
112  C     old_AREA :: hold sea-ice fraction field before any seaice-thermo update  C     old_AREA :: hold sea-ice fraction field before any seaice-thermo update
113        _RL old_AREA     (1:sNx,1:sNy)        _RL old_AREA     (1:sNx,1:sNy)
114  # ifdef SEAICE_AGE_VOL  # ifdef SEAICE_AGE_VOL
115  C     old_HEFF :: hold sea-ice effective thicness field before any seaice-thermo update  C     old_HEFF :: hold sea-ice effective thicness field before any seaice-thermo update
116        _RL old_HEFF     (1:sNx,1:sNy)        _RL old_HEFF     (1:sNx,1:sNy)
117        _RL age_actual        _RL age_actual
118  # endif /* SEAICE_AGE_VOL */  # endif /* SEAICE_AGE_VOL */
119  #endif /* SEAICE_AGE */  #endif /* SEAICE_AGE */
120    
121  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
122        _RL DIAGarray     (1:sNx,1:sNy)        _RL DIAGarray     (1:sNx,1:sNy)
123        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
124        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
125  #endif  #endif
126    
127        IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN        IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
128         kSurface        = Nr         kSurface        = Nr
129        ELSE        ELSE
130         kSurface        = 1         kSurface        = 1
131        ENDIF        ENDIF
132    
133  C     FREEZING TEMP. OF SEA WATER (deg C)  C     FREEZING TEMP. OF SEA WATER (deg C)
134        TBC          = SEAICE_freeze        TBC          = SEAICE_freeze
135  C     RATIO OF WATER DESITY TO SNOW DENSITY  C     RATIO OF WATER DESITY TO SNOW DENSITY
136        SDF          = 1000.0 _d 0/SEAICE_rhoSnow        SDF          = 1000.0 _d 0/SEAICE_rhoSnow
137  C     RATIO OF SEA ICE DENSITY to SNOW DENSITY  C     RATIO OF SEA ICE DENSITY to SNOW DENSITY
138        ICE2SNOW     = SDF * ICE2WATR        ICE2SNOW     = SDF * ICE2WATR
139  C     HEAT OF FUSION OF ICE (J/m^3)  C     HEAT OF FUSION OF ICE (J/m^3)
140        QI           = 302.0 _d +06        QI           = 302.0 _d +06
141        recip_QI     = 1.0 _d 0 / QI        recip_QI     = 1.0 _d 0 / QI
142  C     HEAT OF FUSION OF SNOW (J/m^3)  C     HEAT OF FUSION OF SNOW (J/m^3)
143        QS           = 1.1 _d +08        QS           = 1.1 _d +08
144    
145        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
146         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
147  c  c
148  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
149            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
150            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
151            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
152            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
153            act3 = myThid - 1            act3 = myThid - 1
154            max3 = nTx*nTy            max3 = nTx*nTy
155            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
156            iicekey = (act1 + 1) + act2*max1            iicekey = (act1 + 1) + act2*max1
157       &                      + act3*max1*max2       &                      + act3*max1*max2
158       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
159  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
160  C  C
161  C     initialise a few fields  C     initialise a few fields
162  C  C
163  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
164  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
165  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
166  CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,  CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,
167  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
168  CADJ STORE qsw(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE qsw(:,:,bi,bj)  = comlev1_bibj,
169  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
170  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
171          DO J=1,sNy          DO J=1,sNy
172           DO I=1,sNx           DO I=1,sNx
173            FHEFF(I,J)      = 0.0 _d 0            FHEFF(I,J)      = 0.0 _d 0
174            FICE (I,J)      = 0.0 _d 0            FICE (I,J)      = 0.0 _d 0
175            QNETO(I,J)      = 0.0 _d 0            QNETO(I,J)      = 0.0 _d 0
176            QNETI(I,J)      = 0.0 _d 0            QNETI(I,J)      = 0.0 _d 0
177            QSWO (I,J)      = 0.0 _d 0            QSWO (I,J)      = 0.0 _d 0
178            QSWI (I,J)      = 0.0 _d 0            QSWI (I,J)      = 0.0 _d 0
179            HCORR(I,J)      = 0.0 _d 0            HCORR(I,J)      = 0.0 _d 0
180            RESID_HEAT(I,J) = 0.0 _d 0            RESID_HEAT(I,J) = 0.0 _d 0
181  #ifdef SEAICE_SALINITY  #ifdef SEAICE_SALINITY
182            saltFluxAdjust(I,J) = 0.0 _d 0            saltFluxAdjust(I,J) = 0.0 _d 0
183  #endif  #endif
184  #ifdef SEAICE_MULTICATEGORY  #ifdef SEAICE_MULTICATEGORY
185            FICEP(I,J)      = 0.0 _d 0            FICEP(I,J)      = 0.0 _d 0
186            QSWIP(I,J)      = 0.0 _d 0            QSWIP(I,J)      = 0.0 _d 0
187  #endif  #endif
188           ENDDO           ENDDO
189          ENDDO          ENDDO
190          DO J=1-oLy,sNy+oLy          DO J=1-oLy,sNy+oLy
191           DO I=1-oLx,sNx+oLx           DO I=1-oLx,sNx+oLx
192            saltWtrIce(I,J,bi,bj) = 0.0 _d 0            saltWtrIce(I,J,bi,bj) = 0.0 _d 0
193            frWtrIce(I,J,bi,bj)   = 0.0 _d 0            frWtrIce(I,J,bi,bj)   = 0.0 _d 0
194  #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION  #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
195            frWtrAtm(I,J,bi,bj)   = 0.0 _d 0            frWtrAtm(I,J,bi,bj)   = 0.0 _d 0
196  #endif  #endif
197           ENDDO           ENDDO
198          ENDDO          ENDDO
199    
200  #ifdef SEAICE_AGE  #ifdef SEAICE_AGE
201  C     store the initial ice fraction over the domain  C     store the initial ice fraction over the domain
202          DO J=1,sNy          DO J=1,sNy
203           DO I=1,sNx           DO I=1,sNx
204             old_AREA(i,j) = AREA(I,J,bi,bj)             old_AREA(i,j) = AREA(I,J,bi,bj)
205  # ifdef SEAICE_AGE_VOL  # ifdef SEAICE_AGE_VOL
206             old_HEFF(i,j) = HEFF(I,J,bi,bj)             old_HEFF(i,j) = HEFF(I,J,bi,bj)
207  # endif  # endif
208           ENDDO           ENDDO
209          ENDDO          ENDDO
210  #endif /* SEAICE_AGE */  #endif /* SEAICE_AGE */
211    
212    
213  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
214  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,
215  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
216  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
217  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
218  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
219          DO J=1,sNy          DO J=1,sNy
220           DO I=1,sNx           DO I=1,sNx
221  C     COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM  C     COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM
222  C     ON ICE THICKNESS FOR BUDGET COMPUTATION  C     ON ICE THICKNESS FOR BUDGET COMPUTATION
223  C     The default of A22 = 0.15 is a common threshold for defining  C     The default of A22 = 0.15 is a common threshold for defining
224  C     the ice edge. This ice concentration usually does not occur  C     the ice edge. This ice concentration usually does not occur
225  C     due to thermodynamics but due to advection.  C     due to thermodynamics but due to advection.
226            areaLoc      = MAX(A22,AREANm1(I,J,bi,bj))            areaLoc      = MAX(A22,AREANm1(I,J,bi,bj))
227            HICE(I,J)    = HEFFNm1(I,J,bi,bj)/areaLoc            HICE(I,J)    = HEFFNm1(I,J,bi,bj)/areaLoc
228  C     Do we know what this is for?  C     Do we know what this is for?
229            HICE(I,J)    = MAX(HICE(I,J),0.05 _d +00)            HICE(I,J)    = MAX(HICE(I,J),0.05 _d +00)
230  C     Capping the actual ice thickness effectively enforces a  C     Capping the actual ice thickness effectively enforces a
231  C     minimum of heat flux through the ice and helps getting rid of  C     minimum of heat flux through the ice and helps getting rid of
232  C     very thick ice.  C     very thick ice.
233  cdm actually, this does exactly the opposite, i.e., ice is thicker  cdm actually, this does exactly the opposite, i.e., ice is thicker
234  cdm when HICE is capped, so I am commenting out  cdm when HICE is capped, so I am commenting out
235  cdm          HICE(I,J)    = MIN(HICE(I,J),9.0 _d +00)  cdm          HICE(I,J)    = MIN(HICE(I,J),9.0 _d +00)
236            hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc            hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc
237           ENDDO           ENDDO
238          ENDDO          ENDDO
239    
240  C NOW DETERMINE MIXED LAYER TEMPERATURE  C NOW DETERMINE MIXED LAYER TEMPERATURE
241          DO J=1,sNy          DO J=1,sNy
242           DO I=1,sNx           DO I=1,sNx
243            TMIX(I,J,bi,bj)=theta(I,J,kSurface,bi,bj)+273.16 _d +00            TMIX(I,J,bi,bj)=theta(I,J,kSurface,bi,bj)+273.16 _d +00
244  #ifdef SEAICE_DEBUG  #ifdef SEAICE_DEBUG
245            TMIX(I,J,bi,bj)=MAX(TMIX(I,J,bi,bj),271.2 _d +00)            TMIX(I,J,bi,bj)=MAX(TMIX(I,J,bi,bj),271.2 _d +00)
246  #endif  #endif
247           ENDDO           ENDDO
248          ENDDO          ENDDO
249    
250  C THERMAL WIND OF ATMOSPHERE  C THERMAL WIND OF ATMOSPHERE
251          DO J=1,sNy          DO J=1,sNy
252           DO I=1,sNx           DO I=1,sNx
253  C     copy the wind speed computed in exf_wind.F to UG  C     copy the wind speed computed in exf_wind.F to UG
254            UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))            UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
255  CML   this is the old code, which does the same  CML   this is the old code, which does the same
256  CML          SPEED_SQ = UWIND(I,J,bi,bj)**2 + VWIND(I,J,bi,bj)**2  CML          SPEED_SQ = UWIND(I,J,bi,bj)**2 + VWIND(I,J,bi,bj)**2
257  CML          IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN  CML          IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
258  CML             UG(I,J)=SEAICE_EPS  CML             UG(I,J)=SEAICE_EPS
259  CML          ELSE  CML          ELSE
260  CML             UG(I,J)=SQRT(SPEED_SQ)  CML             UG(I,J)=SQRT(SPEED_SQ)
261  CML          ENDIF  CML          ENDIF
262           ENDDO           ENDDO
263          ENDDO          ENDDO
264    
265    
266  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
267  cphCADJ STORE heff   = comlev1, key = ikey_dynamics, byte = isbyte  cphCADJ STORE heff   = comlev1, key = ikey_dynamics, byte = isbyte
268  cphCADJ STORE hsnow  = comlev1, key = ikey_dynamics, byte = isbyte  cphCADJ STORE hsnow  = comlev1, key = ikey_dynamics, byte = isbyte
269  cphCADJ STORE uwind  = comlev1, key = ikey_dynamics, byte = isbyte  cphCADJ STORE uwind  = comlev1, key = ikey_dynamics, byte = isbyte
270  cphCADJ STORE vwind  = comlev1, key = ikey_dynamics, byte = isbyte  cphCADJ STORE vwind  = comlev1, key = ikey_dynamics, byte = isbyte
271  c  c
272  CADJ STORE tice   = comlev1, key = ikey_dynamics, byte = isbyte  CADJ STORE tice   = comlev1, key = ikey_dynamics, byte = isbyte
273  # ifdef SEAICE_MULTICATEGORY  # ifdef SEAICE_MULTICATEGORY
274  CADJ STORE tices  = comlev1, key = ikey_dynamics, byte = isbyte  CADJ STORE tices  = comlev1, key = ikey_dynamics, byte = isbyte
275  # endif  # endif
276  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
277    
278  C NOW DETERMINE GROWTH RATES  C NOW DETERMINE GROWTH RATES
279  C FIRST DO OPEN WATER  C FIRST DO OPEN WATER
280          CALL SEAICE_BUDGET_OCEAN(          CALL SEAICE_BUDGET_OCEAN(
281       I       UG,       I       UG,
282       U       TMIX,       U       TMIX,
283       O       QNETO, QSWO,       O       QNETO, QSWO,
284       I       bi, bj, myTime, myIter, myThid )       I       bi, bj, myTime, myIter, myThid )
285    
286  C NOW DO ICE  C NOW DO ICE
287          IF (useRelativeWind) THEN          IF (useRelativeWind) THEN
288  C     Compute relative wind speed over sea ice.  C     Compute relative wind speed over sea ice.
289           DO J=1,sNy           DO J=1,sNy
290            DO I=1,sNx            DO I=1,sNx
291             SPEED_SQ =             SPEED_SQ =
292       &          (uWind(I,J,bi,bj)       &          (uWind(I,J,bi,bj)
293       &          +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)       &          +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
294       &                    +uVel(i+1,j,kSurface,bi,bj))       &                    +uVel(i+1,j,kSurface,bi,bj))
295       &          -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2       &          -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
296       &          +(vWind(I,J,bi,bj)       &          +(vWind(I,J,bi,bj)
297       &          +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)       &          +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
298       &                    +vVel(i,j+1,kSurface,bi,bj))       &                    +vVel(i,j+1,kSurface,bi,bj))
299       &          -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2       &          -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
300             IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN             IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
301               UG(I,J)=SEAICE_EPS               UG(I,J)=SEAICE_EPS
302             ELSE             ELSE
303               UG(I,J)=SQRT(SPEED_SQ)               UG(I,J)=SQRT(SPEED_SQ)
304             ENDIF             ENDIF
305            ENDDO            ENDDO
306           ENDDO           ENDDO
307          ENDIF          ENDIF
308  #ifdef SEAICE_MULTICATEGORY  #ifdef SEAICE_MULTICATEGORY
309  C--  Start loop over muli-categories  C--  Start loop over muli-categories
310          DO IT=1,MULTDIM          DO IT=1,MULTDIM
311  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
312           ilockey = (iicekey-1)*MULTDIM + IT           ilockey = (iicekey-1)*MULTDIM + IT
313  CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,  CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,
314  CADJ &                           key = ilockey, byte = isbyte  CADJ &                           key = ilockey, byte = isbyte
315  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
316           RK=REAL(IT)           RK=REAL(IT)
317           DO J=1,sNy           DO J=1,sNy
318            DO I=1,sNx            DO I=1,sNx
319             HICEP(I,J)=(HICE(I,J)/MULTDIM)*((2.0 _d 0*RK)-1.0 _d 0)             HICEP(I,J)=(HICE(I,J)/MULTDIM)*((2.0 _d 0*RK)-1.0 _d 0)
320             TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)             TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)
321            ENDDO            ENDDO
322           ENDDO           ENDDO
323           CALL SEAICE_BUDGET_ICE(           CALL SEAICE_BUDGET_ICE(
324       I        UG, HICEP, hSnwLoc,       I        UG, HICEP, hSnwLoc,
325       U        TICE,       U        TICE,
326       O        FICEP, QSWIP,       O        FICEP, QSWIP,
327       I        bi, bj, myTime, myIter, myThid )       I        bi, bj, myTime, myIter, myThid )
328           DO J=1,sNy           DO J=1,sNy
329            DO I=1,sNx            DO I=1,sNx
330  C     average surface heat fluxes/growth rates  C     average surface heat fluxes/growth rates
331             FICE (I,J)          = FICE(I,J) + FICEP(I,J)/MULTDIM             FICE (I,J)          = FICE(I,J) + FICEP(I,J)/MULTDIM
332             QSWI (I,J)          = QSWI(I,J) + QSWIP(I,J)/MULTDIM             QSWI (I,J)          = QSWI(I,J) + QSWIP(I,J)/MULTDIM
333             TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)             TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)
334            ENDDO            ENDDO
335           ENDDO           ENDDO
336          ENDDO          ENDDO
337  C--  End loop over multi-categories  C--  End loop over multi-categories
338  #else  /* SEAICE_MULTICATEGORY */  #else  /* SEAICE_MULTICATEGORY */
339          CALL SEAICE_BUDGET_ICE(          CALL SEAICE_BUDGET_ICE(
340       I       UG, HICE, hSnwLoc,       I       UG, HICE, hSnwLoc,
341       U       TICE,       U       TICE,
342       O       FICE, QSWI,       O       FICE, QSWI,
343       I       bi, bj, myTime, myIter, myThid )       I       bi, bj, myTime, myIter, myThid )
344  #endif /* SEAICE_MULTICATEGORY */  #endif /* SEAICE_MULTICATEGORY */
345    
346  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
347          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
348           IF ( DIAGNOSTICS_IS_ON('SIatmQnt',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIatmQnt',myThid) ) THEN
349            DO J=1,sNy            DO J=1,sNy
350             DO I=1,sNx             DO I=1,sNx
351              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (
352       &           FICE(I,J) * areaNm1(I,J,bi,bj) +       &           FICE(I,J) * areaNm1(I,J,bi,bj) +
353       &           QNETO(I,J) * ( ONE - areaNm1(I,J,bi,bj) ) )       &           QNETO(I,J) * ( ONE - areaNm1(I,J,bi,bj) ) )
354             ENDDO             ENDDO
355            ENDDO            ENDDO
356            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)
357           ENDIF           ENDIF
358          ENDIF          ENDIF
359  #endif  #endif
360    
361  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
362  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,
363  CADJ &                          key = iicekey, byte = isbyte  CADJ &                          key = iicekey, byte = isbyte
364  CADJ STORE heff(:,:,bi,bj)    = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)    = comlev1_bibj,
365  CADJ &                          key = iicekey, byte = isbyte  CADJ &                          key = iicekey, byte = isbyte
366  #endif  #endif
367  C  C
368  C--   compute and apply ice growth due to oceanic heat flux from below  C--   compute and apply ice growth due to oceanic heat flux from below
369  C  C
370          DO J=1,sNy          DO J=1,sNy
371           DO I=1,sNx           DO I=1,sNx
372  C--   Create or melt sea-ice so that first-level oceanic temperature  C--   Create or melt sea-ice so that first-level oceanic temperature
373  C     is approximately at the freezing point when there is sea-ice.  C     is approximately at the freezing point when there is sea-ice.
374  C     Initially the units of YNEG/availHeat are m of sea-ice.  C     Initially the units of YNEG/availHeat are m of sea-ice.
375  C     The factor dRf(1)/72.0764, used to convert temperature  C     The factor dRf(1)/72.0764, used to convert temperature
376  C     change in deg K to m of sea-ice, is approximately:  C     change in deg K to m of sea-ice, is approximately:
377  C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)  C     dRf(1) * (sea water heat capacity = 3996 J/kg/K)
378  C        * (density of sea-water = 1026 kg/m^3)  C        * (density of sea-water = 1026 kg/m^3)
379  C        / (latent heat of fusion of sea-ice = 334000 J/kg)  C        / (latent heat of fusion of sea-ice = 334000 J/kg)
380  C        / (density of sea-ice = 910 kg/m^3)  C        / (density of sea-ice = 910 kg/m^3)
381  C     Negative YNEG/availHeat leads to ice growth.  C     Negative YNEG/availHeat leads to ice growth.
382  C     Positive YNEG/availHeat leads to ice melting.  C     Positive YNEG/availHeat leads to ice melting.
383            IF ( .NOT. inAdMode ) THEN            IF ( .NOT. inAdMode ) THEN
384  #ifdef SEAICE_VARIABLE_FREEZING_POINT  #ifdef SEAICE_VARIABLE_FREEZING_POINT
385             TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0             TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
386  #endif /* SEAICE_VARIABLE_FREEZING_POINT */  #endif /* SEAICE_VARIABLE_FREEZING_POINT */
387             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN             IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
388                availHeat(i,j) = SEAICE_availHeatFrac                availHeat(i,j) = SEAICE_availHeatFrac
389       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
390       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
391             ELSE             ELSE
392                availHeat(i,j) = SEAICE_availHeatFracFrz                availHeat(i,j) = SEAICE_availHeatFracFrz
393       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)       &             * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
394       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0       &             * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
395             ENDIF             ENDIF
396            ELSE            ELSE
397             availHeat(i,j) = 0.             availHeat(i,j) = 0.
398            ENDIF            ENDIF
399  C     local copy of old effective ice thickness  C     local copy of old effective ice thickness
400            hEffOld(I,J)      = HEFF(I,J,bi,bj)            hEffOld(I,J)      = HEFF(I,J,bi,bj)
401  C     Melt (availHeat>0) or create (availHeat<0) sea ice  C     Melt (availHeat>0) or create (availHeat<0) sea ice
402            HEFF(I,J,bi,bj) = MAX(ZERO,HEFF(I,J,bi,bj)-availHeat(I,J))            HEFF(I,J,bi,bj) = MAX(ZERO,HEFF(I,J,bi,bj)-availHeat(I,J))
403  C  C
404            YNEG(I,J,bi,bj)   = hEffOld(I,J)    - HEFF(I,J,bi,bj)            YNEG(I,J,bi,bj)   = hEffOld(I,J)    - HEFF(I,J,bi,bj)
405  C  C
406            saltWtrIce(I,J,bi,bj)   = saltWtrIce(I,J,bi,bj)            saltWtrIce(I,J,bi,bj)   = saltWtrIce(I,J,bi,bj)
407       &                              - YNEG(I,J,bi,bj)       &                              - YNEG(I,J,bi,bj)
408            RESID_HEAT(I,J)   = availHeat(I,J)  - YNEG(I,J,bi,bj)            RESID_HEAT(I,J)   = availHeat(I,J)  - YNEG(I,J,bi,bj)
409  C     YNEG now contains m of ice melted (>0) or created (<0)  C     YNEG now contains m of ice melted (>0) or created (<0)
410  C     saltWtrIce contains m of ice melted (<0) or created (>0)  C     saltWtrIce contains m of ice melted (<0) or created (>0)
411  C     RESID_HEAT is residual heat above freezing in equivalent m of ice  C     RESID_HEAT is residual heat above freezing in equivalent m of ice
412           ENDDO           ENDDO
413          ENDDO          ENDDO
414    
415  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
416          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
417           IF ( DIAGNOSTICS_IS_ON('SIyneg  ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIyneg  ',myThid) ) THEN
418            CALL DIAGNOSTICS_FILL(YNEG,'SIyneg  ',0,1,1,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(YNEG,'SIyneg  ',0,1,1,bi,bj,myThid)
419           ENDIF           ENDIF
420          ENDIF          ENDIF
421  #endif  #endif
422    
423  cph(  cph(
424  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
425  cphCADJ STORE heff   = comlev1, key = ikey_dynamics, byte = isbyte  cphCADJ STORE heff   = comlev1, key = ikey_dynamics, byte = isbyte
426  cphCADJ STORE hsnow  = comlev1, key = ikey_dynamics, byte = isbyte  cphCADJ STORE hsnow  = comlev1, key = ikey_dynamics, byte = isbyte
427  #endif  #endif
428  cph)  cph)
429  c  c
430  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
431  CADJ STORE area(:,:,bi,bj)   = comlev1_bibj,  CADJ STORE area(:,:,bi,bj)   = comlev1_bibj,
432  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
433  CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj)  = comlev1_bibj,
434  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
435  CADJ STORE fice(:,:)         = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
436  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
437  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
438  cph)  cph)
439  C  C
440  C--   compute and apply ice growth due to atmospheric fluxes from above  C--   compute and apply ice growth due to atmospheric fluxes from above
441  C  C
442          DO J=1,sNy          DO J=1,sNy
443           DO I=1,sNx           DO I=1,sNx
444  C NOW CALCULATE CORRECTED effective growth in J/m^2 (>0=melt)  C NOW CALCULATE CORRECTED effective growth in J/m^2 (>0=melt)
445            GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J)*AREANm1(I,J,bi,bj)            GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J)*AREANm1(I,J,bi,bj)
446           ENDDO           ENDDO
447          ENDDO          ENDDO
448    
449  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
450  CADJ STORE fice(:,:)         = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
451  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
452  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
453    
454          DO J=1,sNy          DO J=1,sNy
455           DO I=1,sNx           DO I=1,sNx
456            IF(FICE(I,J).LT.ZERO.AND.AREANm1(I,J,bi,bj).GT.ZERO) THEN            IF(FICE(I,J).LT.ZERO.AND.AREANm1(I,J,bi,bj).GT.ZERO) THEN
457  C     use FICE to melt snow and CALCULATE CORRECTED GROWTH  C     use FICE to melt snow and CALCULATE CORRECTED GROWTH
458  C     effective snow thickness in J/m^2  C     effective snow thickness in J/m^2
459             snowEnergy=HSNOW(I,J,bi,bj)*QS             snowEnergy=HSNOW(I,J,bi,bj)*QS
460             IF(GHEFF(I,J).LE.snowEnergy) THEN             IF(GHEFF(I,J).LE.snowEnergy) THEN
461  C     not enough heat to melt all snow; use up all heat flux FICE  C     not enough heat to melt all snow; use up all heat flux FICE
462              HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS              HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS
463  C     SNOW CONVERTED INTO WATER AND THEN INTO equivalent m of ICE melt  C     SNOW CONVERTED INTO WATER AND THEN INTO equivalent m of ICE melt
464  C     The factor 1/ICE2SNOW converts m of snow to m of sea-ice  C     The factor 1/ICE2SNOW converts m of snow to m of sea-ice
465              frWtrIce(I,J,bi,bj) =              frWtrIce(I,J,bi,bj) =
466       &       frWtrIce(I,J,bi,bj) - GHEFF(I,J)/(QS*ICE2SNOW)       &       frWtrIce(I,J,bi,bj) - GHEFF(I,J)/(QS*ICE2SNOW)
467              FICE    (I,J) = ZERO              FICE    (I,J) = ZERO
468             ELSE             ELSE
469  C     enought heat to melt snow completely;  C     enought heat to melt snow completely;
470  C     compute remaining heat flux that will melt ice  C     compute remaining heat flux that will melt ice
471              FICE(I,J)=-(GHEFF(I,J)-snowEnergy)/              FICE(I,J)=-(GHEFF(I,J)-snowEnergy)/
472       &           SEAICE_deltaTtherm/AREANm1(I,J,bi,bj)       &           SEAICE_deltaTtherm/AREANm1(I,J,bi,bj)
473  C     convert all snow to melt water (fresh water flux)  C     convert all snow to melt water (fresh water flux)
474              frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)              frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)
475       &           -HSNOW(I,J,bi,bj)/ICE2SNOW       &           -HSNOW(I,J,bi,bj)/ICE2SNOW
476              HSNOW(I,J,bi,bj)=0.0 _d 0              HSNOW(I,J,bi,bj)=0.0 _d 0
477             END IF             END IF
478            END IF            END IF
479           ENDDO           ENDDO
480          ENDDO          ENDDO
481    
482  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
483  CADJ STORE fice(:,:)         = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
484  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
485  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
486    
487          DO J=1,sNy          DO J=1,sNy
488           DO I=1,sNx           DO I=1,sNx
489  C     now get cell averaged growth rate in W/m^2, >0 causes ice growth  C     now get cell averaged growth rate in W/m^2, >0 causes ice growth
490            FHEFF(I,J)= FICE(I,J)  * AREANm1(I,J,bi,bj)            FHEFF(I,J)= FICE(I,J)  * AREANm1(I,J,bi,bj)
491       &              + QNETO(I,J) * (ONE-AREANm1(I,J,bi,bj))       &              + QNETO(I,J) * (ONE-AREANm1(I,J,bi,bj))
492           ENDDO           ENDDO
493          ENDDO          ENDDO
494  cph(  cph(
495  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
496  CADJ STORE heff(:,:,bi,bj)   = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj)   = comlev1_bibj,
497  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
498  CADJ STORE fice(:,:)         = comlev1_bibj,  CADJ STORE fice(:,:)         = comlev1_bibj,
499  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
500  CADJ STORE fheff(:,:)        = comlev1_bibj,  CADJ STORE fheff(:,:)        = comlev1_bibj,
501  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
502  CADJ STORE qneto(:,:)        = comlev1_bibj,  CADJ STORE qneto(:,:)        = comlev1_bibj,
503  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
504  CADJ STORE qswi(:,:)         = comlev1_bibj,  CADJ STORE qswi(:,:)         = comlev1_bibj,
505  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
506  CADJ STORE qswo(:,:)         = comlev1_bibj,  CADJ STORE qswo(:,:)         = comlev1_bibj,
507  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
508  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
509  cph)  cph)
510  C  C
511  C     First update (freeze or melt) ice area  C     First update (freeze or melt) ice area
512  C  C
513          DO J=1,sNy          DO J=1,sNy
514           DO I=1,sNx           DO I=1,sNx
515  C     negative growth in meters of ice (>0 for melting)  C     negative growth in meters of ice (>0 for melting)
516            growthNeg  = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI            growthNeg  = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI
517  C     negative growth must not exceed effective ice thickness (=volume)  C     negative growth must not exceed effective ice thickness (=volume)
518  C     (that is, cannot melt more than all the ice)  C     (that is, cannot melt more than all the ice)
519            growthHEFF = -ONE*MIN(HEFF(I,J,bi,bj),growthNeg)            growthHEFF = -ONE*MIN(HEFF(I,J,bi,bj),growthNeg)
520  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
521            DIAGarray(I,J) = growthHEFF            DIAGarray(I,J) = growthHEFF
522  #endif  #endif
523  C     growthHEFF < 0 means melting  C     growthHEFF < 0 means melting
524            HCORR(I,J) = MIN(ZERO,growthHEFF)            HCORR(I,J) = MIN(ZERO,growthHEFF)
525  C     gain of new effective ice thickness over open water (>0 by definition)  C     gain of new effective ice thickness over open water (>0 by definition)
526            GAREA(I,J) = MAX(ZERO,SEAICE_deltaTtherm*QNETO(I,J)*recip_QI)            GAREA(I,J) = MAX(ZERO,SEAICE_deltaTtherm*QNETO(I,J)*recip_QI)
527  CML   removed these loops and moved TAMC store directive up  CML   removed these loops and moved TAMC store directive up
528  CML         ENDDO  CML         ENDDO
529  CML        ENDDO  CML        ENDDO
530  CML#ifdef ALLOW_AUTODIFF_TAMC  CML#ifdef ALLOW_AUTODIFF_TAMC
531  CMLCADJ STORE area(:,:,bi,bj) = comlev1_bibj,  CMLCADJ STORE area(:,:,bi,bj) = comlev1_bibj,
532  CMLCADJ &                       key = iicekey, byte = isbyte  CMLCADJ &                       key = iicekey, byte = isbyte
533  CML#endif  CML#endif
534  CML        DO J=1,sNy  CML        DO J=1,sNy
535  CML         DO I=1,sNx  CML         DO I=1,sNx
536  C     Here we finally compute the new AREA  C     Here we finally compute the new AREA
537            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
538             AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+             AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+
539       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO_south       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO_south
540       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)
541       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)
542            ELSE            ELSE
543             AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+             AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+
544       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO       &          (ONE-AREANm1(I,J,bi,bj))*GAREA(I,J)/HO
545       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)       &          +HALF*HCORR(I,J)*AREANm1(I,J,bi,bj)
546       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)       &          /(HEFF(I,J,bi,bj)+.00001 _d 0)
547            ENDIF            ENDIF
548           ENDDO           ENDDO
549          ENDDO          ENDDO
550    
551  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
552          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
553           IF ( DIAGNOSTICS_IS_ON('SIfice  ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIfice  ',myThid) ) THEN
554            CALL DIAGNOSTICS_FILL(DIAGarray,'SIfice  ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIfice  ',0,1,3,bi,bj,myThid)
555           ENDIF           ENDIF
556          ENDIF          ENDIF
557  #endif  #endif
558    
559  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
560  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
561  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
562  #endif  #endif
563  C  C
564  C     now update (freeze or melt) HEFF  C     now update (freeze or melt) HEFF
565  C  C
566          DO J=1,sNy          DO J=1,sNy
567           DO I=1,sNx           DO I=1,sNx
568  C     negative growth (>0 for melting) of existing ice in meters  C     negative growth (>0 for melting) of existing ice in meters
569            growthNeg       = -SEAICE_deltaTtherm*            growthNeg       = -SEAICE_deltaTtherm*
570       &                      FICE(I,J)*recip_QI*AREANm1(I,J,bi,bj)       &                      FICE(I,J)*recip_QI*AREANm1(I,J,bi,bj)
571  C     negative growth must not exceed effective ice thickness (=volume)  C     negative growth must not exceed effective ice thickness (=volume)
572  C     (that is, cannot melt more than all the ice)  C     (that is, cannot melt more than all the ice)
573            growthHEFF      = -ONE*MIN(HEFF(I,J,bi,bj),growthNeg)            growthHEFF      = -ONE*MIN(HEFF(I,J,bi,bj),growthNeg)
574  C     growthHEFF < 0 means melting  C     growthHEFF < 0 means melting
575            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + growthHEFF            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + growthHEFF
576  C     add effective growth to fresh water of ice  C     add effective growth to fresh water of ice
577            saltWtrIce(I,J,bi,bj) = saltWtrIce(I,J,bi,bj) + growthHEFF            saltWtrIce(I,J,bi,bj) = saltWtrIce(I,J,bi,bj) + growthHEFF
578    
579  C     now calculate QNETI under ice (if any) as the difference between  C     now calculate QNETI under ice (if any) as the difference between
580  C     the available "heat flux" growthNeg and the actual growthHEFF;  C     the available "heat flux" growthNeg and the actual growthHEFF;
581  C     keep in mind that growthNeg and growthHEFF have different signs  C     keep in mind that growthNeg and growthHEFF have different signs
582  C     by construction  C     by construction
583            QNETI(I,J) =  (growthHEFF + growthNeg)*QI/SEAICE_deltaTtherm            QNETI(I,J) =  (growthHEFF + growthNeg)*QI/SEAICE_deltaTtherm
584    
585  C     now update other things  C     now update other things
586    
587  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
588            IF(FICE(I,J).GT.ZERO) THEN            IF(FICE(I,J).GT.ZERO) THEN
589  C     freezing, add precip as snow  C     freezing, add precip as snow
590             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*
591       &            PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)*SDF       &            PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)*SDF
592            ELSE            ELSE
593  C     add precip as rain, water converted into equivalent m of  C     add precip as rain, water converted into equivalent m of
594  C     ice by 1/ICE2WATR.  C     ice by 1/ICE2WATR.
595  C     Do not get confused by the sign:  C     Do not get confused by the sign:
596  C     precip   > 0 for downward flux of fresh water  C     precip   > 0 for downward flux of fresh water
597  C     frWtrIce > 0 for more ice (corresponds to an upward "fresh water flux"),  C     frWtrIce > 0 for more ice (corresponds to an upward "fresh water flux"),
598  C     so that here the rain is added *as if* it is melted ice (which is not  C     so that here the rain is added *as if* it is melted ice (which is not
599  C     true, but just a trick; physically the rain just runs as water  C     true, but just a trick; physically the rain just runs as water
600  C     through the ice into the ocean)  C     through the ice into the ocean)
601             frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)             frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)
602       &            -PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)*       &            -PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)*
603       &            SEAICE_deltaTtherm/ICE2WATR       &            SEAICE_deltaTtherm/ICE2WATR
604            ENDIF            ENDIF
605  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
606    
607           ENDDO           ENDDO
608          ENDDO          ENDDO
609    
610  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
611  cphCADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  cphCADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
612  cphCADJ &                        key = iicekey, byte = isbyte  cphCADJ &                        key = iicekey, byte = isbyte
613  #endif  #endif
614    
615  cph( very sensitive bit here by JZ  cph( very sensitive bit here by JZ
616  #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING  #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
617          DO J=1,sNy          DO J=1,sNy
618           DO I=1,sNx           DO I=1,sNx
619  C     Now melt snow if there is residual heat left in surface level  C     Now melt snow if there is residual heat left in surface level
620  C     Note that units of YNEG and frWtrIce are m of ice  C     Note that units of YNEG and frWtrIce are m of ice
621            IF( RESID_HEAT(I,J) .GT. ZERO .AND.            IF( RESID_HEAT(I,J) .GT. ZERO .AND.
622       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
623             GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR,             GHEFF(I,J)       = MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR,
624       &                             RESID_HEAT(I,J) )       &                             RESID_HEAT(I,J) )
625             YNEG(I,J,bi,bj)  = YNEG(I,J,bi,bj) +GHEFF(I,J)             YNEG(I,J,bi,bj)  = YNEG(I,J,bi,bj) +GHEFF(I,J)
626            ENDIF            ENDIF
627           ENDDO           ENDDO
628          ENDDO          ENDDO
629    
630  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
631  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
632  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
633  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
634  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
635  #endif  #endif
636          DO J=1,sNy          DO J=1,sNy
637           DO I=1,sNx           DO I=1,sNx
638            IF( RESID_HEAT(I,J) .GT. ZERO .AND.            IF( RESID_HEAT(I,J) .GT. ZERO .AND.
639       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN       &       HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
640             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-GHEFF(I,J)*SDF*ICE2WATR             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-GHEFF(I,J)*SDF*ICE2WATR
641             frWtrIce(I,J,bi,bj)    = frWtrIce(I,J,bi,bj)-GHEFF(I,J)             frWtrIce(I,J,bi,bj)    = frWtrIce(I,J,bi,bj)-GHEFF(I,J)
642            ENDIF            ENDIF
643           ENDDO           ENDDO
644          ENDDO          ENDDO
645  #endif  #endif
646  cph)  cph)
647    
648  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
649  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
650  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
651  # ifdef SEAICE_SALINITY  # ifdef SEAICE_SALINITY
652  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
653  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
654  # endif /* SEAICE_SALINITY */  # endif /* SEAICE_SALINITY */
655  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
656    
657  #ifdef ALLOW_ATM_TEMP  #ifdef ALLOW_ATM_TEMP
658          DO J=1,sNy          DO J=1,sNy
659           DO I=1,sNx           DO I=1,sNx
660  C NOW GET FRESH WATER FLUX  C NOW GET FRESH WATER FLUX
661            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(            EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(
662       &         ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )       &         ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
663       &         * ( ONE - AREANm1(I,J,bi,bj) )       &         * ( ONE - AREANm1(I,J,bi,bj) )
664  #ifdef ALLOW_RUNOFF  #ifdef ALLOW_RUNOFF
665       &         - RUNOFF(I,J,bi,bj)       &         - RUNOFF(I,J,bi,bj)
666  #endif  #endif
667       &         + frWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm       &         + frWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm
668       &         + saltWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm       &         + saltWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm
669       &         )*rhoConstFresh       &         )*rhoConstFresh
670           ENDDO           ENDDO
671          ENDDO          ENDDO
672    
673  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
674          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
675           IF ( DIAGNOSTICS_IS_ON('SIatmFW ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIatmFW ',myThid) ) THEN
676            DO J=1,sNy            DO J=1,sNy
677             DO I=1,sNx             DO I=1,sNx
678              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)*(              DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)*(
679       &           PRECIP(I,J,bi,bj)       &           PRECIP(I,J,bi,bj)
680       &           - EVAP(I,J,bi,bj)       &           - EVAP(I,J,bi,bj)
681       &           *( ONE - AREANm1(I,J,bi,bj) )       &           *( ONE - AREANm1(I,J,bi,bj) )
682       &           + RUNOFF(I,J,bi,bj)       &           + RUNOFF(I,J,bi,bj)
683       &           )*rhoConstFresh       &           )*rhoConstFresh
684             ENDDO             ENDDO
685            ENDDO            ENDDO
686            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmFW ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmFW ',0,1,3,bi,bj,myThid)
687           ENDIF           ENDIF
688          ENDIF          ENDIF
689  #endif  #endif
690  #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION  #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
691          DO J=1,sNy          DO J=1,sNy
692           DO I=1,sNx           DO I=1,sNx
693            frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(            frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
694       &         PRECIP(I,J,bi,bj)       &         PRECIP(I,J,bi,bj)
695       &         - EVAP(I,J,bi,bj)       &         - EVAP(I,J,bi,bj)
696       &         *( ONE - AREANm1(I,J,bi,bj) )       &         *( ONE - AREANm1(I,J,bi,bj) )
697       &         + RUNOFF(I,J,bi,bj)       &         + RUNOFF(I,J,bi,bj)
698       &         )*rhoConstFresh       &         )*rhoConstFresh
699           ENDDO           ENDDO
700          ENDDO          ENDDO
701  #endif  #endif
702    
703  C COMPUTE SURFACE SALT FLUX AND ADJUST ICE SALINITY  C COMPUTE SURFACE SALT FLUX AND ADJUST ICE SALINITY
704    
705  #ifdef SEAICE_SALINITY  #ifdef SEAICE_SALINITY
706    
707          DO J=1,sNy          DO J=1,sNy
708           DO I=1,sNx           DO I=1,sNx
709  C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean  C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
710            IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN            IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
711               saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *               saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
712       &            HSALT(I,J,bi,bj) / SEAICE_deltaTtherm       &            HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
713               HSALT(I,J,bi,bj) = 0.0 _d 0               HSALT(I,J,bi,bj) = 0.0 _d 0
714            ENDIF            ENDIF
715           ENDDO           ENDDO
716          ENDDO          ENDDO
717    
718  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
719  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
720  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
721  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
722    
723          DO J=1,sNy          DO J=1,sNy
724           DO I=1,sNx           DO I=1,sNx
725  C saltWtrIce > 0 : m of sea ice that is created  C saltWtrIce > 0 : m of sea ice that is created
726            IF ( saltWtrIce(I,J,bi,bj) .GE. 0.0 ) THEN            IF ( saltWtrIce(I,J,bi,bj) .GE. 0.0 ) THEN
727               saltFlux(I,J,bi,bj) =               saltFlux(I,J,bi,bj) =
728       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
729       &            ICE2WATR*rhoConstFresh*SEAICE_salinity*       &            ICE2WATR*rhoConstFresh*SEAICE_salinity*
730       &            salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm       &            salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm
731  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
732  C saltPlumeFlux is defined only during freezing:  C saltPlumeFlux is defined only during freezing:
733               saltPlumeFlux(I,J,bi,bj)=               saltPlumeFlux(I,J,bi,bj)=
734       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
735       &            ICE2WATR*rhoConstFresh*(1-SEAICE_salinity)*       &            ICE2WATR*rhoConstFresh*(1-SEAICE_salinity)*
736       &            salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm       &            salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm
737  C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean  C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
738               IF ( .NOT. SaltPlumeSouthernOcean ) THEN               IF ( .NOT. SaltPlumeSouthernOcean ) THEN
739                IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )                IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
740       &             saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0       &             saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
741               ENDIF               ENDIF
742    
743  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
744  C saltWtrIce < 0 : m of sea ice that is melted  C saltWtrIce < 0 : m of sea ice that is melted
745            ELSE            ELSE
746               saltFlux(I,J,bi,bj) =               saltFlux(I,J,bi,bj) =
747       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*       &            HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
748       &            HSALT(I,J,bi,bj)/       &            HSALT(I,J,bi,bj)/
749       &            (HEFF(I,J,bi,bj)-saltWtrIce(I,J,bi,bj))/       &            (HEFF(I,J,bi,bj)-saltWtrIce(I,J,bi,bj))/
750       &            SEAICE_deltaTtherm       &            SEAICE_deltaTtherm
751  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
752               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
753  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
754            ENDIF            ENDIF
755  C update HSALT based on surface saltFlux  C update HSALT based on surface saltFlux
756            HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +            HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
757       &         saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm       &         saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
758            saltFlux(I,J,bi,bj) =            saltFlux(I,J,bi,bj) =
759       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
760  C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean  C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
761            IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN            IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
762               saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -               saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
763       &            HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /       &            HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /
764       &            SEAICE_deltaTtherm       &            SEAICE_deltaTtherm
765               HSALT(I,J,bi,bj) = 0.0 _d 0               HSALT(I,J,bi,bj) = 0.0 _d 0
766  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
767               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0               saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
768  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
769            ENDIF            ENDIF
770           ENDDO           ENDDO
771          ENDDO          ENDDO
772  #endif /* SEAICE_SALINITY */  #endif /* SEAICE_SALINITY */
773  #endif /* ALLOW_ATM_TEMP */  #endif /* ALLOW_ATM_TEMP */
774    
775          DO J=1,sNy          DO J=1,sNy
776           DO I=1,sNx           DO I=1,sNx
777  C NOW GET TOTAL QNET AND QSW  C NOW GET TOTAL QNET AND QSW
778            QNET(I,J,bi,bj) = QNETI(I,J) * AREANm1(I,J,bi,bj)            QNET(I,J,bi,bj) = QNETI(I,J) * AREANm1(I,J,bi,bj)
779       &                     +QNETO(I,J) * (ONE-AREANm1(I,J,bi,bj))       &                     +QNETO(I,J) * (ONE-AREANm1(I,J,bi,bj))
780            QSW(I,J,bi,bj)  = QSWI(I,J)  * AREANm1(I,J,bi,bj)            QSW(I,J,bi,bj)  = QSWI(I,J)  * AREANm1(I,J,bi,bj)
781       &                     +QSWO(I,J)  * (ONE-AREANm1(I,J,bi,bj))       &                     +QSWO(I,J)  * (ONE-AREANm1(I,J,bi,bj))
782           ENDDO           ENDDO
783          ENDDO          ENDDO
784    
785  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
786          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
787           IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN
788            DO J=1,sNy            DO J=1,sNy
789             DO I=1,sNx             DO I=1,sNx
790              DIAGarray(I,J) = QNETO(I,J) * (ONE-areaNm1(I,J,bi,bj))              DIAGarray(I,J) = QNETO(I,J) * (ONE-areaNm1(I,J,bi,bj))
791             ENDDO             ENDDO
792            ENDDO            ENDDO
793            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)
794           ENDIF           ENDIF
795           IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN
796            DO J=1,sNy            DO J=1,sNy
797             DO I=1,sNx             DO I=1,sNx
798              DIAGarray(I,J) = QNETI(I,J) * areaNm1(I,J,bi,bj)              DIAGarray(I,J) = QNETI(I,J) * areaNm1(I,J,bi,bj)
799             ENDDO             ENDDO
800            ENDDO            ENDDO
801            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)
802           ENDIF           ENDIF
803          ENDIF          ENDIF
804  #endif  #endif
805    
806  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
807  cphCADJ STORE yneg(:,:,bi,bj) = comlev1_bibj,  cphCADJ STORE yneg(:,:,bi,bj) = comlev1_bibj,
808  cphCADJ &                       key = iicekey, byte = isbyte  cphCADJ &                       key = iicekey, byte = isbyte
809  #endif  #endif
810          DO J=1,sNy          DO J=1,sNy
811           DO I=1,sNx           DO I=1,sNx
812  C Now convert YNEG back to deg K.  C Now convert YNEG back to deg K.
813            YNEG(I,J,bi,bj) = YNEG(I,J,bi,bj)*recip_dRf(kSurface) *            YNEG(I,J,bi,bj) = YNEG(I,J,bi,bj)*recip_dRf(kSurface) *
814       &                      recip_hFacC(i,j,kSurface,bi,bj)*72.0764 _d 0       &                      recip_hFacC(i,j,kSurface,bi,bj)*72.0764 _d 0
815           ENDDO           ENDDO
816          ENDDO          ENDDO
817    
818  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
819  CADJ STORE yneg(:,:,bi,bj) = comlev1_bibj,  CADJ STORE yneg(:,:,bi,bj) = comlev1_bibj,
820  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
821  #endif  #endif
822          DO J=1,sNy          DO J=1,sNy
823           DO I=1,sNx           DO I=1,sNx
824  C Add YNEG contribution to QNET  C Add YNEG contribution to QNET
825            QNET(I,J,bi,bj) = QNET(I,J,bi,bj)            QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
826       &         +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm       &         +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm
827       &         *maskC(I,J,kSurface,bi,bj)       &         *maskC(I,J,kSurface,bi,bj)
828       &         *HeatCapacity_Cp*rUnit2mass       &         *HeatCapacity_Cp*rUnit2mass
829       &         *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)       &         *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)
830           ENDDO           ENDDO
831          ENDDO          ENDDO
832    
833  #ifdef SEAICE_DEBUG  #ifdef SEAICE_DEBUG
834         CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )         CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
835         CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )         CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
836         CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )         CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
837  #endif /* SEAICE_DEBUG */  #endif /* SEAICE_DEBUG */
838    
839  crg Added by Ralf Giering: do we need DO_WE_NEED_THIS ?  crg Added by Ralf Giering: do we need DO_WE_NEED_THIS ?
840  #define DO_WE_NEED_THIS  #define DO_WE_NEED_THIS
841  C NOW ZERO OUTSIDE POINTS  C NOW ZERO OUTSIDE POINTS
842  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
843  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
844  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
845  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
846  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
847  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
848          DO J=1,sNy          DO J=1,sNy
849           DO I=1,sNx           DO I=1,sNx
850  C NOW SET AREA(I,J,bi,bj)=0 WHERE NO ICE IS  C NOW SET AREA(I,J,bi,bj)=0 WHERE NO ICE IS
851            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj)            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj)
852       &                         ,HEFF(I,J,bi,bj)/.0001 _d 0)       &                         ,HEFF(I,J,bi,bj)/.0001 _d 0)
853           ENDDO           ENDDO
854          ENDDO          ENDDO
855  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
856  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,  CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
857  CADJ &                       key = iicekey, byte = isbyte  CADJ &                       key = iicekey, byte = isbyte
858  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
859          DO J=1,sNy          DO J=1,sNy
860           DO I=1,sNx           DO I=1,sNx
861  C NOW TRUNCATE AREA  C NOW TRUNCATE AREA
862  #ifdef DO_WE_NEED_THIS  #ifdef DO_WE_NEED_THIS
863            AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))            AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
864           ENDDO           ENDDO
865          ENDDO          ENDDO
866  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
867  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
868  CADJ &                        key = iicekey, byte = isbyte  CADJ &                        key = iicekey, byte = isbyte
869  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
870  CADJ &                         key = iicekey, byte = isbyte  CADJ &                         key = iicekey, byte = isbyte
871  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
872          DO J=1,sNy          DO J=1,sNy
873           DO I=1,sNx           DO I=1,sNx
874            AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))            AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
875            HSNOW(I,J,bi,bj)  = MAX(ZERO,HSNOW(I,J,bi,bj))            HSNOW(I,J,bi,bj)  = MAX(ZERO,HSNOW(I,J,bi,bj))
876  #endif /* DO_WE_NEED_THIS */  #endif /* DO_WE_NEED_THIS */
877            AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)            AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
878            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
879  #ifdef SEAICE_CAP_HEFF  #ifdef SEAICE_CAP_HEFF
880            HEFF(I,J,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,bi,bj))            HEFF(I,J,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,bi,bj))
881  #endif /* SEAICE_CAP_HEFF */  #endif /* SEAICE_CAP_HEFF */
882            HSNOW(I,J,bi,bj)  = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)            HSNOW(I,J,bi,bj)  = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
883           ENDDO           ENDDO
884          ENDDO          ENDDO
885    
886  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
887          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
888           IF ( DIAGNOSTICS_IS_ON('SIthdgrh',myThid) ) THEN           IF ( DIAGNOSTICS_IS_ON('SIthdgrh',myThid) ) THEN
889  C     use (abuse) GHEFF to diagnose the total thermodynamic growth rate  C     use (abuse) GHEFF to diagnose the total thermodynamic growth rate
890            DO J=1,sNy            DO J=1,sNy
891             DO I=1,sNx             DO I=1,sNx
892              GHEFF(I,J) = (HEFF(I,J,bi,bj)-HEFFNm1(I,J,bi,bj))              GHEFF(I,J) = (HEFF(I,J,bi,bj)-HEFFNm1(I,J,bi,bj))
893       &           /SEAICE_deltaTtherm       &           /SEAICE_deltaTtherm
894             ENDDO             ENDDO
895            ENDDO            ENDDO
896            CALL DIAGNOSTICS_FILL(GHEFF,'SIthdgrh',0,1,3,bi,bj,myThid)            CALL DIAGNOSTICS_FILL(GHEFF,'SIthdgrh',0,1,3,bi,bj,myThid)
897           ENDIF           ENDIF
898          ENDIF          ENDIF
899  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
900    
901  #ifdef ALLOW_SEAICE_FLOODING  #ifdef ALLOW_SEAICE_FLOODING
902          IF ( SEAICEuseFlooding ) THEN          IF ( SEAICEuseFlooding ) THEN
903  C     convert snow to ice if submerged  C     convert snow to ice if submerged
904           DO J=1,sNy           DO J=1,sNy
905            DO I=1,sNx            DO I=1,sNx
906             hDraft     = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow             hDraft     = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
907       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/1000. _d 0       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/1000. _d 0
908  C     here GEFF is the gain of ice due to flooding  C     here GEFF is the gain of ice due to flooding
909             GHEFF(I,J) = hDraft - MIN(hDraft,HEFF(I,J,bi,bj))             GHEFF(I,J) = hDraft - MIN(hDraft,HEFF(I,J,bi,bj))
910             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + GHEFF(I,J)             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + GHEFF(I,J)
911             HSNOW(I,J,bi,bj)  = MAX(0. _d 0,             HSNOW(I,J,bi,bj)  = MAX(0. _d 0,
912       &          HSNOW(I,J,bi,bj)-GHEFF(I,J)*ICE2SNOW)       &          HSNOW(I,J,bi,bj)-GHEFF(I,J)*ICE2SNOW)
913            ENDDO            ENDDO
914           ENDDO           ENDDO
915  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
916           IF ( useDiagnostics ) THEN           IF ( useDiagnostics ) THEN
917            IF ( DIAGNOSTICS_IS_ON('SIsnwice',myThid) ) THEN            IF ( DIAGNOSTICS_IS_ON('SIsnwice',myThid) ) THEN
918  C     turn GHEFF into a rate  C     turn GHEFF into a rate
919             DO J=1,sNy             DO J=1,sNy
920              DO I=1,sNx              DO I=1,sNx
921               GHEFF(I,J) = GHEFF(I,J)/SEAICE_deltaTtherm               GHEFF(I,J) = GHEFF(I,J)/SEAICE_deltaTtherm
922              ENDDO              ENDDO
923             ENDDO             ENDDO
924             CALL DIAGNOSTICS_FILL(GHEFF,'SIsnwice',0,1,3,bi,bj,myThid)             CALL DIAGNOSTICS_FILL(GHEFF,'SIsnwice',0,1,3,bi,bj,myThid)
925            ENDIF            ENDIF
926           ENDIF           ENDIF
927  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
928          ENDIF          ENDIF
929  #endif /* ALLOW_SEAICE_FLOODING */  #endif /* ALLOW_SEAICE_FLOODING */
930    
931          IF ( useRealFreshWaterFlux ) THEN          IF ( useRealFreshWaterFlux ) THEN
932           DO J=1,sNy           DO J=1,sNy
933            DO I=1,sNx            DO I=1,sNx
934             sIceLoad(i,j,bi,bj) = HEFF(I,J,bi,bj)*SEAICE_rhoIce             sIceLoad(i,j,bi,bj) = HEFF(I,J,bi,bj)*SEAICE_rhoIce
935       &                         + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow       &                         + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
936            ENDDO            ENDDO
937           ENDDO           ENDDO
938          ENDIF          ENDIF
939    
940  #ifdef SEAICE_AGE  #ifdef SEAICE_AGE
941  # ifndef SEAICE_AGE_VOL  # ifndef SEAICE_AGE_VOL
942  C     Sources and sinks for sea ice age:  C     Sources and sinks for sea ice age:
943  C     assume that a) freezing: new ice fraction forms with zero age  C     assume that a) freezing: new ice fraction forms with zero age
944  C                 b) melting: ice fraction vanishes with current age  C                 b) melting: ice fraction vanishes with current age
945          DO J=1,sNy          DO J=1,sNy
946           DO I=1,sNx           DO I=1,sNx
947            IF ( AREA(I,J,bi,bj) .GT. 0.15 ) THEN            IF ( AREA(I,J,bi,bj) .GT. 0.15 ) THEN
948             IF ( AREA(i,j,bi,bj) .LT. old_AREA(i,j) ) THEN             IF ( AREA(i,j,bi,bj) .LT. old_AREA(i,j) ) THEN
949  C--   scale effective ice-age to account for ice-age sink associated with melting  C--   scale effective ice-age to account for ice-age sink associated with melting
950              IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)              IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
951       &         *AREA(i,j,bi,bj)/old_AREA(i,j)       &         *AREA(i,j,bi,bj)/old_AREA(i,j)
952             ENDIF             ENDIF
953  C--   account for aging:  C--   account for aging:
954             IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)             IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
955       &        + AREA(i,j,bi,bj) * SEAICE_deltaTtherm       &        + AREA(i,j,bi,bj) * SEAICE_deltaTtherm
956            ELSE            ELSE
957             IceAge(i,j,bi,bj) = ZERO             IceAge(i,j,bi,bj) = ZERO
958            ENDIF            ENDIF
959           ENDDO           ENDDO
960          ENDDO          ENDDO
961  # else /* ifdef SEAICE_AGE_VOL */  # else /* ifdef SEAICE_AGE_VOL */
962  C     Sources and sinks for sea ice age:  C     Sources and sinks for sea ice age:
963  C     assume that a) freezing: new ice volume forms with zero age  C     assume that a) freezing: new ice volume forms with zero age
964  C                 b) melting: ice volume vanishes with current age  C                 b) melting: ice volume vanishes with current age
965          DO J=1,sNy          DO J=1,sNy
966           DO I=1,sNx           DO I=1,sNx
967  C--   compute actual age from effective age:  C--   compute actual age from effective age:
968            IF (OLD_AREA(i,j).GT.0. _d 0) THEN            IF (OLD_AREA(i,j).GT.0. _d 0) THEN
969             age_actual=IceAge(i,j,bi,bj)/OLD_AREA(i,j)             age_actual=IceAge(i,j,bi,bj)/OLD_AREA(i,j)
970            ELSE            ELSE
971             age_actual=0. _d 0             age_actual=0. _d 0
972            ENDIF                    ENDIF        
973            IF ( (OLD_HEFF(i,j).LT.HEFF(i,j,bi,bj)).AND.            IF ( (OLD_HEFF(i,j).LT.HEFF(i,j,bi,bj)).AND.
974       &         (AREA(i,j,bi,bj).GT.0.15) ) THEN       &         (AREA(i,j,bi,bj).GT.0.15) ) THEN
975             age_actual=age_actual*OLD_HEFF(i,j)/             age_actual=age_actual*OLD_HEFF(i,j)/
976       &          HEFF(i,j,bi,bj)+SEAICE_deltaTtherm       &          HEFF(i,j,bi,bj)+SEAICE_deltaTtherm
977            ELSEIF (AREA(i,j,bi,bj).LE.0.15) THEN            ELSEIF (AREA(i,j,bi,bj).LE.0.15) THEN
978             age_actual=0. _d 0             age_actual=0. _d 0
979            ELSE            ELSE
980             age_actual=age_actual+SEAICE_deltaTtherm             age_actual=age_actual+SEAICE_deltaTtherm
981            ENDIF            ENDIF
982  C--   re-scale to effective age:  C--   re-scale to effective age:
983            IceAge(i,j,bi,bj) = age_actual*AREA(i,j,bi,bj)            IceAge(i,j,bi,bj) = age_actual*AREA(i,j,bi,bj)
984           ENDDO           ENDDO
985          ENDDO          ENDDO
986  # endif /* SEAICE_AGE_VOL */  # endif /* SEAICE_AGE_VOL */
987  #endif /* SEAICE_AGE */  #endif /* SEAICE_AGE */
988    
989         ENDDO         ENDDO
990        ENDDO        ENDDO
991    
992        RETURN        RETURN
993        END        END
Legend:
Removed from v.1.68  
changed lines
  Added in v.1.69
  ViewVC Help
Powered by ViewVC 1.1.22