/[MITgcm]/MITgcm/pkg/thsice/thsice_extend.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_extend.F

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

revision 1.3 by jmc, Fri Feb 10 00:30:32 2006 UTC revision 1.4 by jmc, Thu May 25 18:03:24 2006 UTC
# Line 7  CBOP Line 7  CBOP
7  C     !ROUTINE: THSICE_EXTEND  C     !ROUTINE: THSICE_EXTEND
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE THSICE_EXTEND(        SUBROUTINE THSICE_EXTEND(
10       I                  esurp, Tf,       I                  bi, bj, siLo, siHi, sjLo, sjHi,
11       U                  sst, compact, iceThick, snowThick, qicen,       I                  iMin,iMax, jMin,jMax, dBugFlag,
12       O                  qleft, fresh, fsalt,       I                  fzMlOc, tFrz, tOce,
13       I                  dBugFlag, myThid )       U                  icFrac, hIce, hSnow,
14         U                  tSrf, tIc1, tIc2, qIc1, qIc2,
15         O                  flx2oc, frw2oc, fsalt,
16         I                  myTime, myIter, myThid )
17  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
18  C     *==========================================================*  C     *==========================================================*
19  C     | S/R  THSICE_EXTEND                                                C     | S/R  THSICE_EXTEND
20  C     | o Extend sea-ice area incresing ice fraction  C     | o Extend sea-ice area incresing ice fraction
21  C     *==========================================================*  C     *==========================================================*
22  C     | o incorporate surplus of energy to  C     | o incorporate surplus of energy to
23  C     |   make new ice or make ice grow laterally  C     |   make new ice or make ice grow laterally
24  C     *==========================================================*  C     *==========================================================*
25  C     \ev  C     \ev
26    
# Line 31  C     == Global variables == Line 34  C     == Global variables ==
34    
35  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
36  C     == Routine Arguments ==  C     == Routine Arguments ==
37  C     esurp    :: energy available for freezing           [W/m2]  C     siLo,siHi   :: size of input/output array: 1rst dim. lower,higher bounds
38  C     Tf       :: freezing temperature  [oC]  C     sjLo,sjHi   :: size of input/output array: 2nd  dim. lower,higher bounds
39  C     sst      :: Sea Surf Temp. [oC]  C     bi,bj       :: tile indices
40  C     compact  :: fraction of grid area covered in ice  C     iMin,iMax   :: computation domain: 1rst index range
41  C     iceThick :: ice height  [m]  C     jMin,jMax   :: computation domain: 2nd  index range
42  C     snowThick:: snow height [m]  C     dBugFlag    :: allow to print debugging stuff (e.g. on 1 grid point).
43  C     qicen    :: ice enthalpy [J/kg]  C---  Input:
44  C     qleft    :: (additional) heat flux to ocean         [W/m2]  C         iceMask :: sea-ice fractional mask [0-1]
45  C     fsalt    :: (additional) salt flux to ocean        [ g/m2/s]  C  fzMlOc (esurp) :: ocean mixed-layer freezing/melting potential [W/m2]
46  C     fresh    :: (additional) fresh water flux to ocean [kg/m2/s]  C  tFrz    (Tf)   :: sea-water freezing temperature [oC] (function of S)
47  C     dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).  C  tOce    (sst)  :: surface level oceanic temperature [oC]
48  C     myThid   :: thread number for this instance of the routine.  C---  Modified (input&output):
49    C  icFrac(iceFrac):: fraction of grid area covered in ice
50    C  hIce (iceThick):: ice height [m]
51    C hSnow(snowThick):: snow height [m]
52    C  tSrf           :: surface (ice or snow) temperature [oC]
53    C  tIc1           :: temperature of ice layer 1 [oC]
54    C  tIc2           :: temperature of ice layer 2 [oC]
55    C  qIc1   (qicen) :: ice enthalpy (J/kg), 1rst level
56    C  qIc2   (qicen) :: ice enthalpy (J/kg), 2nd level
57    C---  Output
58    C  flx2oc   (=)   :: (additional) heat flux to ocean    [W/m2]        (+=dwn)
59    C  frw2oc   (=)   :: (additional) fresh water flux to ocean [kg/m2/s] (+=dwn)
60    C  fsalt    (=)   :: (additional) salt flux to ocean        [g/m2/s]  (+=dwn)
61    C---  Input:
62    C     myTime      :: current Time of simulation [s]
63    C     myIter      :: current Iteration number in simulation
64    C     myThid      :: my Thread Id number
65          INTEGER siLo, siHi, sjLo, sjHi
66          INTEGER bi,bj
67          INTEGER iMin, iMax
68          INTEGER jMin, jMax
69          LOGICAL dBugFlag
70    c     _RL iceMask(siLo:siHi,sjLo:sjHi)
71          _RL fzMlOc (siLo:siHi,sjLo:sjHi)
72          _RL tFrz   (siLo:siHi,sjLo:sjHi)
73          _RL tOce   (siLo:siHi,sjLo:sjHi)
74          _RL icFrac (siLo:siHi,sjLo:sjHi)
75          _RL hIce   (siLo:siHi,sjLo:sjHi)
76          _RL hSnow  (siLo:siHi,sjLo:sjHi)
77          _RL tSrf   (siLo:siHi,sjLo:sjHi)
78          _RL tIc1   (siLo:siHi,sjLo:sjHi)
79          _RL tIc2   (siLo:siHi,sjLo:sjHi)
80          _RL qIc1   (siLo:siHi,sjLo:sjHi)
81          _RL qIc2   (siLo:siHi,sjLo:sjHi)
82          _RL flx2oc (siLo:siHi,sjLo:sjHi)
83          _RL frw2oc (siLo:siHi,sjLo:sjHi)
84          _RL fsalt  (siLo:siHi,sjLo:sjHi)
85          _RL  myTime
86          INTEGER myIter
87          INTEGER myThid
88    CEOP
89    
90    #ifdef ALLOW_THSICE
91    C     !LOCAL VARIABLES:
92    C---  local copy of input/output argument list variables (see description above)
93        _RL esurp        _RL esurp
94        _RL Tf        _RL Tf
95        _RL sst        _RL sst
96        _RL compact        _RL iceFrac
97        _RL iceThick        _RL iceThick
98        _RL snowThick        _RL snowThick
99        _RL qicen(nlyr)        _RL qicen(nlyr)
       _RL qleft  
       _RL fresh  
       _RL fsalt  
       LOGICAL dBugFlag  
       INTEGER myThid  
 CEOP  
100    
 #ifdef ALLOW_THSICE  
   
 C     !LOCAL VARIABLES:  
101  C     == Local variables ==  C     == Local variables ==
102  C     qicAv    :: mean enthalpy of ice (layer 1 & 2) [J/m^3]  C     qicAv    :: mean enthalpy of ice (layer 1 & 2) [J/m^3]
103        _RL deltaTice ! time-step for ice model        _RL deltaTice ! time-step for ice model
104        _RL newIce        _RL newIce
105        _RL newIceFrac        _RL newIceFrac
       _RL iceFraction  
106        _RL qicAv        _RL qicAv
107        LOGICAL dBug        INTEGER  i,j     ! loop indices
108    
109    C-    define grid-point location where to print debugging values
110    #include "THSICE_DEBUG.h"
111    
112   1010 FORMAT(A,I3,3F8.3)   1010 FORMAT(A,I3,3F8.3)
113   1020 FORMAT(A,1P4E11.3)   1020 FORMAT(A,1P4E11.3)
       dBug = .FALSE.  
 c     dBug = dBugFlag  
114    
115    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
116          deltaTice = thSIce_deltaT
117    
118    
119          DO j = jMin, jMax
120           DO i = iMin, iMax
121            IF (fzMlOc(i,j).GT.0. _d 0) THEN
122              esurp   = fzMlOc(i,j)
123              Tf      = tFrz(i,j)
124              sst     = tOce(i,j)
125              iceFrac = icFrac(i,j)
126              iceThick= hIce(i,j)
127              snowThick=hSnow(i,j)
128              qicen(1)= qIc1(i,j)
129              qicen(2)= qIc2(i,j)
130    C---
131  C--   start ice  C--   start ice
         deltaTice = thSIce_deltaT  
         iceFraction = compact  
132          newIceFrac = 0. _d 0          newIceFrac = 0. _d 0
133    
134  C-    enthalpy of new ice to form :  C-    enthalpy of new ice to form :
135          IF ( compact.LE.0. _d 0 ) THEN          IF ( iceFrac.LE.0. _d 0 ) THEN
136            qicen(1)= -cpwater*Tmlt1            qicen(1)= -cpwater*Tmlt1
137       &             + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)       &             + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)
138            qicen(2)= -cpice *Tf + Lfresh            qicen(2)= -cpice *Tf + Lfresh
# Line 88  C-    enthalpy of new ice to form : Line 140  C-    enthalpy of new ice to form :
140          qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0          qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0
141          newIce = esurp*deltaTice/qicAv          newIce = esurp*deltaTice/qicAv
142    
143          IF (iceFraction.EQ.0. _d 0) THEN          IF (icFrac(i,j).EQ.0. _d 0) THEN
144  c         IF (newIce.GE.himin*iceMaskmax) THEN  c         IF (newIce.GE.himin*iceMaskmax) THEN
145  C- jmc: above is the original version, but below seems more logical:  C- jmc: above is the original version, but below seems more logical:
146            IF (newIce.GE.himin0*iceMaskmin) THEN            IF (newIce.GE.himin0*iceMaskmin) THEN
147  C-    if there is no ice in grid and enough ice to form:  C-    if there is no ice in grid and enough ice to form:
148              iceThick   = MAX(himin0,newIce/iceMaskmax)              iceThick   = MAX(himin0,newIce/iceMaskmax)
149              newIceFrac = MIN(newIce/himin0,iceMaskmax)              newIceFrac = MIN(newIce/himin0,iceMaskmax)
150              compact = newIceFrac              iceFrac = newIceFrac
151              sst=Tf              sst=Tf
152            ENDIF            ENDIF
153          ELSE          ELSE
154  C-    if there is already some ice  C-    if there is already some ice
155              newIceFrac=MIN(newIce/iceThick,iceMaskmax-iceFraction)              newIceFrac=MIN(newIce/iceThick,iceMaskmax-icFrac(i,j))
156              compact = iceFraction + newIceFrac              iceFrac = icFrac(i,j) + newIceFrac
157  C-    spread snow out over ice  C-    spread snow out over ice
158              snowThick = snowThick*iceFraction/compact              snowThick = snowThick*icFrac(i,j)/iceFrac
159              sst=(1. _d 0-newIceFrac)*sst+newIceFrac*Tf                    sst=(1. _d 0-newIceFrac)*sst+newIceFrac*Tf
160          ENDIF          ENDIF
161          qleft= iceThick*newIceFrac*qicAv/deltaTice  C-    oceanic fluxes:
162          fresh=-(rhoi*iceThick)*newIceFrac/deltaTice          flx2oc(i,j)= iceThick*newIceFrac*qicAv/deltaTice
163          fsalt=-(rhoi*iceThick*saltice)*newIceFrac/deltaTice          frw2oc(i,j)=-(rhoi*iceThick)*newIceFrac/deltaTice
164            fsalt(i,j)= -(rhoi*iceThick*saltice)*newIceFrac/deltaTice
165          IF (dBug) WRITE(6,1020) 'ThSI_EXT: iceH, newIce, newIceFrac=',  
166       &     iceThick, newIce, newIceFrac  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
167    #ifdef ALLOW_DBUG_THSICE
168            IF ( dBug(i,j,bi,bj) ) THEN
169              WRITE(6,1020) 'ThSI_EXT: iceH, newIce, newIceFrac=',
170         &                             iceThick, newIce, newIceFrac
171              WRITE(6,1020) 'ThSI_EXT: iceFrac,flx2oc,fsalt,frw2oc=',
172         &                  iceFrac,flx2oc(i,j),fsalt(i,j),frw2oc(i,j)
173            ENDIF
174    #endif
175    #ifdef CHECK_ENERGY_CONSERV
176              CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 1,
177         I            icFrac(i,j), iceFrac, iceThick, snowThick, qicen,
178         I            flx2oc(i,j), frw2oc(i,j), fsalt(i,j),
179         I            myTime, myIter, myThid )
180    #endif /* CHECK_ENERGY_CONSERV */
181    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
182    C-- Update Sea-Ice state output:
183              IF ( iceFrac.GT.0. _d 0 .AND. icFrac(i,j).EQ.0. _d 0) THEN
184                tSrf(i,j)   = tFrz(i,j)
185                tIc1(i,j)   = tFrz(i,j)
186                tIc2(i,j)   = tFrz(i,j)
187                qIc1(i,j)   = qicen(1)
188                qIc2(i,j)   = qicen(2)
189              ENDIF
190              icFrac(i,j) = iceFrac
191              hIce(i,j)   = iceThick
192              hSnow(i,j ) = snowThick
193            ENDIF
194           ENDDO
195          ENDDO
196    
197  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
198    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22