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

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

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

revision 1.3 by jmc, Wed Apr 7 23:40:34 2004 UTC revision 1.30 by gforget, Fri Dec 17 04:00:14 2010 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "THSICE_OPTIONS.h"  #include "THSICE_OPTIONS.h"
5    #ifdef ALLOW_ATM2D
6    # include "ctrparam.h"
7    #endif
8    
9  CBOP  CBOP
10  C     !ROUTINE: THSICE_STEP_FWD  C     !ROUTINE: THSICE_STEP_FWD
11  C     !INTERFACE:  C     !INTERFACE:
12        SUBROUTINE THSICE_STEP_FWD(        SUBROUTINE THSICE_STEP_FWD(
13       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
14       I             prcAtm,       I             prcAtm,
      U             evpAtm, flxSW,  
15       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | S/R  THSICE_STEP_FWD              C     | S/R  THSICE_STEP_FWD
19  C     | o Step Forward Therm-SeaIce model.  C     | o Step Forward Therm-SeaIce model.
20  C     *==========================================================*  C     *==========================================================*
21  C     \ev  C     \ev
# Line 26  C     === Global variables === Line 28  C     === Global variables ===
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "FFIELDS.h"  #include "FFIELDS.h"
31    #ifdef  ALLOW_ATM2D
32    # include "ATMSIZE.h"
33    # include "ATM2D_VARS.h"
34    #endif
35  #include "THSICE_SIZE.h"  #include "THSICE_SIZE.h"
36  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
37  #include "THSICE_VARS.h"  #include "THSICE_VARS.h"
38  #include "THSICE_TAVE.h"  #include "THSICE_TAVE.h"
39    #ifdef ALLOW_AUTODIFF_TAMC
40    # include "tamc.h"
41    # include "tamc_keys.h"
42    #endif
43    
44          INTEGER siLo, siHi, sjLo, sjHi
45          PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
46          PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
47    
48  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
49  C     === Routine arguments ===  C     === Routine arguments ===
50    C- input:
51  C     bi,bj   :: tile indices  C     bi,bj   :: tile indices
52  C   iMin,iMax :: computation domain: 1rst index range  C   iMin,iMax :: computation domain: 1rst index range
53  C   jMin,jMax :: computation domain: 2nd  index range  C   jMin,jMax :: computation domain: 2nd  index range
54    C     prcAtm  :: total precip from the atmosphere [kg/m2/s]
55    C     myTime  :: current Time of simulation [s]
56    C     myIter  :: current Iteration number in simulation
57    C     myThid  :: my Thread Id number
58    C-- Use fluxes hold in commom blocks
59  C- input:  C- input:
60  C     prcAtm  :: total precip from the atmosphere [kg/m2/s]  C     icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
61  C     evpAtm  :: (Inp) evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)  C     icFlxAtm  :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
62  C     flxSW   :: (Inp) short-wave heat flux (+=down): downward comp. only  C     icFrwAtm  :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
 C                      (part.1), becomes net SW flux into ocean (part.2).  
63  C- output  C- output
64  C     evpAtm  :: (Out) net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)  C     icFlxAtm  :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
65  C     flxSW   :: (Out) net surf. heat flux from the atmosphere [W/m2], (+=down)  C     icFrwAtm  :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
 C     myTime  :: time counter for this thread  
 C     myIter  :: iteration counter for this thread  
 C     myThid  :: thread number for this instance of the routine.  
66        INTEGER bi,bj        INTEGER bi,bj
67        INTEGER iMin, iMax        INTEGER iMin, iMax
68        INTEGER jMin, jMax        INTEGER jMin, jMax
69        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
70        _RL  myTime        _RL  myTime
71        INTEGER myIter        INTEGER myIter
72        INTEGER myThid        INTEGER myThid
# Line 61  CEOP Line 75  CEOP
75  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
76  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
77  C     === Local variables ===  C     === Local variables ===
78  C     snowPr    :: snow precipitation [kg/m2/s]  C     iceFrac   :: fraction of grid area covered in ice
 C     agingTime :: aging time scale (s)  
 C     ageFac    :: snow aging factor [1]  
 C     albedo    :: surface albedo [0-1]  
 C     flxAtm    :: net heat flux from the atmosphere (+=down) [W/m2]  
 C     frwAtm    :: net fresh-water flux (E-P) to the atmosphere  [kg/m2/s]  
 C     Fbot      :: the oceanic heat flux already incorporated (ice_therm)  
79  C     flx2oc    :: net heat flux from the ice to the ocean (+=down) [W/m2]  C     flx2oc    :: net heat flux from the ice to the ocean (+=down) [W/m2]
80  C     frw2oc    :: fresh-water flux from the ice to the ocean  C     frw2oc    :: fresh-water flux from the ice to the ocean
81  C     fsalt     :: mass salt flux to the ocean  C     fsalt     :: mass salt flux to the ocean
82  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
83  C     TFrzOce   :: sea-water freezing temperature [oC] (function of S)  C     tFrzOce   :: sea-water freezing temperature [oC] (function of S)
84        INTEGER i,j  C     isIceFree :: true for ice-free grid-cell that remains ice-free
85        _RL snowPr  C     ageFac    :: snow aging factor [1]
86        _RL agingTime, ageFac  C     snowFac   :: snowing refreshing-age factor [units of 1/snowPr]
87        _RL albedo        LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     iceFrac  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _RL frwAtm        _RL     flx2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RL flx2oc        _RL     frw2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91        _RL frw2oc        _RL     fsalt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL fsalt        _RL     tFrzOce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93        _RL TFrzOce, cphm, frzmltMxL        _RL     frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL Fbot, esurp        _RL ageFac
95          _RL snowFac
96          _RL cphm
97        _RL opFrac, icFrac        _RL opFrac, icFrac
98        _RL oceV2s, oceTs  #ifdef ALLOW_DIAGNOSTICS
99        _RL compact, hIce, hSnow, Tsf, Tice(nlyr), qicen(nlyr)        _RL tmpFac
100        _RL tmpflx(0:2), tmpdTs  #endif
101          INTEGER i,j
102        LOGICAL dBug        LOGICAL dBugFlag
103    
104        dBug = .FALSE.  C-    define grid-point location where to print debugging values
105   1010 FORMAT(A,1P4E11.3)  #include "THSICE_DEBUG.h"
106    
107        IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN   1010 FORMAT(A,1P4E14.6)
        DO j = jMin, jMax  
         DO i = iMin, iMax  
 c        dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.15 )  
108    
109  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
 C    part.1 : ice-covered fraction ;  
 C     Solve for surface and ice temperature (implicitly) ; compute surf. fluxes  
 C-------  
          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN  
           icFrac  = iceMask(i,j,bi,bj)  
           TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)  
           hIce    = iceHeight(i,j,bi,bj)  
           hSnow   = snowHeight(i,j,bi,bj)  
           Tsf     = Tsrf(i,j,bi,bj)  
           qicen(1)= Qice1(i,j,bi,bj)  
           qicen(2)= Qice2(i,j,bi,bj)  
   
           CALL THSICE_ALBEDO(  
      I               hIce, hSnow, Tsf, snowAge(i,j,bi,bj),  
      O               albedo,  
      I               myThid )  
           flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo)  
   
           CALL THSICE_SOLVE4TEMP(  
      I          useBulkforce, tmpflx, TFrzOce, hIce, hSnow,  
      U          flxSW(i,j), Tsf, qicen,  
      O          Tice, sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),  
      O          tmpdTs, flxAtm(i,j), evpAtm(i,j),  
      I          i,j, bi,bj, myThid)  
110    
111  #ifdef SHORTWAVE_HEATING  #ifdef ALLOW_AUTODIFF_TAMC
112  C--    Update Fluxes :        act1 = bi - myBxLo(myThid)
113            opFrac= 1. _d 0-icFrac        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
114            Qsw(i,j,bi,bj)=-icFrac*flxSW(i,j) +opFrac*Qsw(i,j,bi,bj)        act2 = bj - myByLo(myThid)
115          max2 = myByHi(myThid) - myByLo(myThid) + 1
116          act3 = myThid - 1
117          max3 = nTx*nTy
118          act4 = ikey_dynamics - 1
119          ticekey = (act1 + 1) + act2*max1
120         &                     + act3*max1*max2
121         &                     + act4*max1*max2*max3
122    #endif /* ALLOW_AUTODIFF_TAMC */
123    
124    C-    Initialise
125          dBugFlag = debugLevel.GE.debLevB
126          DO j = 1-OLy, sNy+OLy
127            DO i = 1-OLx, sNx+OLx
128              isIceFree(i,j) = .FALSE.
129    #ifdef ALLOW_ATM2D
130              sFluxFromIce(i,j) = 0. _d 0
131    #else
132              saltFlux(i,j,bi,bj) = 0. _d 0
133    #endif
134    #ifdef ALLOW_AUTODIFF_TAMC
135              iceFrac(i,j) = 0.
136  #endif  #endif
 C--    Update Sea-Ice state :  
           Tsrf(i,j,bi,bj) =Tsf  
           Tice1(i,j,bi,bj)=Tice(1)  
           Tice2(i,j,bi,bj)=Tice(2)  
           Qice1(i,j,bi,bj)=qicen(1)  
           Qice2(i,j,bi,bj)=qicen(2)  
 #ifdef ALLOW_TIMEAVE  
           ice_albedo_Ave(i,j,bi,bj) = ice_albedo_Ave(i,j,bi,bj)  
      &                              + icFrac*albedo*thSIce_deltaT  
 #endif /*ALLOW_TIMEAVE*/  
          ENDIF  
137          ENDDO          ENDDO
138          ENDDO
139    
140          ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
141          snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
142    
143    #ifdef ALLOW_AUTODIFF_TAMC
144    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
145    CADJ STORE iceheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
146    CADJ STORE icfrwatm(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
147    CADJ STORE qice1(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
148    CADJ STORE qice2(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
149    CADJ STORE snowheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
150    #endif
151          DO j = jMin, jMax
152           DO i = iMin, iMax
153            IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
154    C--  Snow aging :
155              snowAge(i,j,bi,bj) = thSIce_deltaT
156         &                       + snowAge(i,j,bi,bj)*ageFac
157              IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
158         &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
159         &          * EXP( - snowFac*snowPrc(i,j,bi,bj) )
160    c    &          * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
161    c    &                  /hNewSnowAge )
162    C-------
163    C note: Any flux of mass (here fresh water) that enter or leave the system
164    C       with a non zero energy HAS TO be counted: add snow precip.
165              icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
166         &                        - Lfresh*snowPrc(i,j,bi,bj)
167    C--
168            ENDIF
169         ENDDO         ENDDO
170          ENDDO
171    
172    #ifdef ALLOW_DIAGNOSTICS
173          IF ( useDiagnostics ) THEN
174            tmpFac = 1. _d 0
175            CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
176            CALL DIAGNOSTICS_FRACT_FILL(
177         I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',
178         I                   0,1,1,bi,bj,myThid)
179            CALL DIAGNOSTICS_FRACT_FILL(
180         I                   siceAlb,   iceMask,tmpFac,1,'SIalbedo',
181         I                   0,1,1,bi,bj,myThid)
182        ENDIF        ENDIF
183        dBug = .FALSE.  #endif /* ALLOW_DIAGNOSTICS */
184          DO j = jMin, jMax
185           DO i = iMin, iMax
186              siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
187           ENDDO
188          ENDDO
189    
190  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
191  C    part.2 : ice-covered fraction ;  C    part.2 : ice-covered fraction ;
192  C     change in ice/snow thickness and ice-fraction  C     change in ice/snow thickness and ice-fraction
193  C     note: can only reduce the ice-fraction but not increase it.  C     note: can only reduce the ice-fraction but not increase it.
194  C-------  C-------
       agingTime = 50. _d 0 * 86400. _d 0  
       ageFac = 1. _d 0 - thSIce_deltaT/agingTime  
195        DO j = jMin, jMax        DO j = jMin, jMax
196         DO i = iMin, iMax         DO i = iMin, iMax
 c       dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.15 )  
197    
198          TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)          tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
         oceTs   = tOceMxL(i,j,bi,bj)  
199          cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)          cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)
200          frzmltMxL = (TFrzOce-oceTs)*cphm/ocean_deltaT          frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
201         &                 * cphm/ocean_deltaT
202          Fbot   = 0. _d 0          iceFrac(i,j) = iceMask(i,j,bi,bj)
203          saltFlux(i,j,bi,bj) = 0. _d 0          flx2oc(i,j)  = icFlxSW(i,j,bi,bj)
         compact= iceMask(i,j,bi,bj)  
204  C-------  C-------
205          IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact.GT.0.) ) THEN  #ifdef ALLOW_DBUG_THSICE
206            WRITE(6,1010) 'ThSI_FWD:-1- iceMask,hIc,hSn,Qnet=',          IF ( dBug(i,j,bi,bj) ) THEN
207       &                  compact, hIce, hSnow, Qnet(i,j,bi,bj)           IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
208            WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL=',            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
209       &                            oceTs,TFrzOce,frzmltMxL            WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
210         &                  iceFrac(i,j), iceHeight(i,j,bi,bj),
211         &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
212              WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
213         &                     tOceMxL(i,j,bi,bj), tFrzOce(i,j),
214         &                     frzmltMxL(i,j), Qnet(i,j,bi,bj)
215             ENDIF
216             IF (iceFrac(i,j).GT.0.)
217         &    WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
218         &      iceFrac(i,j), icFlxAtm(i,j,bi,bj),
219         &      icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
220          ENDIF          ENDIF
221  C-------  #endif
222          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN         ENDDO
223          ENDDO
           oceV2s  = v2ocMxL(i,j,bi,bj)  
           snowPr  = snowPrc(i,j,bi,bj)  
           hIce    = iceHeight(i,j,bi,bj)  
           hSnow   = snowHeight(i,j,bi,bj)  
           Tsf     = Tsrf(i,j,bi,bj)  
           qicen(1)= Qice1(i,j,bi,bj)  
           qicen(2)= Qice2(i,j,bi,bj)  
           flx2oc  = flxSW(i,j)  
   
           CALL THSICE_CALC_THICKN(  
      I          frzmltMxL, TFrzOce, oceTs, oceV2s, snowPr,  
      I          sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj), evpAtm(i,j),  
      U          compact, hIce, hSnow, Tsf, qicen, flx2oc,  
      O          frw2oc, fsalt, Fbot,  
      I          dBug, myThid)  
   
 C- note : snowPr was not supposed to be modified in THSICE_THERM ;  
 C         but to reproduce old results, is reset to zero if Tsf >= 0  
           snowPrc(i,j,bi,bj) = snowPr  
   
 C--  Snow aging :  
           snowAge(i,j,bi,bj) = thSIce_deltaT  
      &                       + snowAge(i,j,bi,bj)*ageFac  
           IF ( snowPr.GT.0. _d 0 )  
      &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)  
      &          * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge )  
 C--  
   
 C-- Diagnostic of Atmospheric Fluxes over sea-ice :  
           frwAtm = evpAtm(i,j) - prcAtm(i,j)  
 C note: Any flux of mass (here fresh water) that enter or leave the system  
 C       with a non zero energy HAS TO be counted: add snow precip.  
           flxAtm(i,j) = flxAtm(i,j) - Lfresh*snowPrc(i,j,bi,bj)  
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
       IF (dBug) WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',  
      &  iceMask(i,j,bi,bj),flxAtm(i,j),evpAtm(i,j),-Lfresh*snowPr  
       IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc=',  
      &   compact,flx2oc,fsalt,frw2oc  
 #ifdef CHECK_ENERGY_CONSERV  
           icFrac = iceMask(i,j,bi,bj)  
           CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,  
      I            icFrac, compact, hIce, hSnow, qicen,  
      I            flx2oc, frw2oc, fsalt, flxAtm(i,j), frwAtm,  
      I            myTime, myIter, myThid )  
 #endif /* CHECK_ENERGY_CONSERV */  
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
224    
225  C--    Update Sea-Ice state :  #ifdef ALLOW_AUTODIFF_TAMC
226  c         iceMask(i,j,bi,bj)=compact  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
227            iceheight(i,j,bi,bj) = hIce  #endif
           snowheight(i,j,bi,bj)= hSnow  
           Tsrf(i,j,bi,bj) =Tsf  
           Qice1(i,j,bi,bj)=qicen(1)  
           Qice2(i,j,bi,bj)=qicen(2)  
228    
229          CALL THSICE_CALC_THICKN(
230         I          bi, bj,
231         I          iMin,iMax, jMin,jMax, dBugFlag,
232         I          iceMask(siLo,sjLo,bi,bj), tFrzOce,
233         I          tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
234         I          snowPrc(siLo,sjLo,bi,bj), prcAtm,
235         I          sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
236         U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
237         U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
238         U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
239         U          icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
240         O          frw2oc, fsalt,
241         I          myTime, myIter, myThid )
242    
243    #ifdef ALLOW_AUTODIFF_TAMC
244    CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
245    CADJ STORE fsalt(:,:)  = comlev1_bibj,key=ticekey,byte=isbyte
246    CADJ STORE flx2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
247    CADJ STORE frw2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
248    #endif
249  C--    Net fluxes :  C--    Net fluxes :
250            frw2oc = frw2oc + (prcAtm(i,j)-snowPrc(i,j,bi,bj))        DO j = jMin, jMax
251           DO i = iMin, iMax
252    c#ifdef ALLOW_AUTODIFF_TAMC
253    c         ikey_1 = i
254    c    &         + sNx*(j-1)
255    c    &         + sNx*sNy*act1
256    c    &         + sNx*sNy*max1*act2
257    c    &         + sNx*sNy*max1*max2*act3
258    c    &         + sNx*sNy*max1*max2*max3*act4
259    c#endif /* ALLOW_AUTODIFF_TAMC */
260    c#ifdef ALLOW_AUTODIFF_TAMC
261    cCADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
262    c#endif
263            IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
264  C-     weighted average net fluxes:  C-     weighted average net fluxes:
265    c#ifdef ALLOW_AUTODIFF_TAMC
266    cCADJ STORE  fsalt(i,j) = comlev1_thsice_1, key=ikey_1
267    cCADJ STORE  flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
268    cCADJ STORE  frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
269    cCADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
270    c#endif
271            icFrac = iceMask(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
272            opFrac= 1. _d 0-icFrac            opFrac= 1. _d 0-icFrac
273            flxAtm(i,j) = icFrac*flxAtm(i,j) - opFrac*Qnet(i,j,bi,bj)  #ifdef ALLOW_ATM2D
274            frwAtm =     icFrac*frwAtm + opFrac*rhofw*EmPmR(i,j,bi,bj)            pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
275            Qnet(i,j,bi,bj)=-icFrac*flx2oc +opFrac*Qnet(i,j,bi,bj)            pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
276            EmPmR(i,j,bi,bj)=-icFrac*frw2oc/rhofw+opFrac*EmPmR(i,j,bi,bj)            sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
277            saltFlux(i,j,bi,bj)=-icFrac*fsalt  #else
278              icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
279            IF (dBug) WRITE(6,1010)'ThSI_FWD:-3- compact,hIc,hSn,Qnet=',       &                        - opFrac*Qnet(i,j,bi,bj)
280       &                      compact,hIce,hSnow,Qnet(i,j,bi,bj)            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
281         &                        + opFrac*EmPmR(i,j,bi,bj)
282          ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN            Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
283            flxAtm(i,j) =  -Qnet(i,j,bi,bj)            EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)
284            frwAtm = rhofw*EmPmR(i,j,bi,bj)       &                    +  opFrac*EmPmR(i,j,bi,bj)
285              saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
286    #endif
287    
288    #ifdef ALLOW_DBUG_THSICE
289              IF (dBug(i,j,bi,bj)) WRITE(6,1010)
290         &          'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
291         &           iceFrac(i,j), iceHeight(i,j,bi,bj),
292         &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
293    #endif
294    
295            ELSEIF (hOceMxL(i,j,bi,bj).GT.0. _d 0) THEN
296              icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
297              icFrwAtm(i,j,bi,bj) = EmPmR(i,j,bi,bj)
298          ELSE          ELSE
299            flxAtm(i,j) = 0. _d 0            icFlxAtm(i,j,bi,bj) = 0. _d 0
300            frwAtm      = 0. _d 0            icFrwAtm(i,j,bi,bj) = 0. _d 0
301          ENDIF          ENDIF
302           ENDDO
303          ENDDO
304    
305  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
306  C    part.3 : freezing of sea-water  C    part.3 : freezing of sea-water
307  C     over ice-free fraction and what is left from ice-covered fraction  C     over ice-free fraction and what is left from ice-covered fraction
308  C-------  C-------
309  c       compact= iceMask(i,j,bi,bj)        CALL THSICE_EXTEND(
310          hIce   = iceHeight(i,j,bi,bj)       I          bi, bj,
311          hSnow  = snowHeight(i,j,bi,bj)       I          iMin,iMax, jMin,jMax, dBugFlag,
312         I          frzmltMxL, tFrzOce,
313          esurp  = frzmltMxL - Fbot*iceMask(i,j,bi,bj)       I          tOceMxL(siLo,sjLo,bi,bj),
314          IF (esurp.GT.0. _d 0) THEN       U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
315            icFrac = compact       U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
316            qicen(1)= Qice1(i,j,bi,bj)       U          Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
317            qicen(2)= Qice2(i,j,bi,bj)       U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
318            CALL THSICE_EXTEND(       O          flx2oc, frw2oc, fsalt,
319       I               esurp, TFrzOce,       I          myTime, myIter, myThid )
320       U               oceTs, compact, hIce, hSnow, qicen,  
321       O               flx2oc, frw2oc, fsalt,  #ifdef ALLOW_AUTODIFF_TAMC
322       I               dBug, myThid )  CADJ STORE snowHeight(:,:,bi,bj) =
323  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  CADJ &     comlev1_bibj,key=ticekey,byte=isbyte
324        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc='  #endif
325       &                 ,compact,flx2oc,fsalt,frw2oc        DO j = jMin, jMax
326  #ifdef CHECK_ENERGY_CONSERV         DO i = iMin, iMax
327            tmpflx(1) = 0.          IF (frzmltMxL(i,j).GT.0. _d 0) THEN
           tmpflx(2) = 0.  
           CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,  
      I            icFrac, compact, hIce, hSnow, qicen,  
      I            flx2oc, frw2oc, fsalt, tmpflx(1), tmpflx(2),  
      I            myTime, myIter, myThid )  
 #endif /* CHECK_ENERGY_CONSERV */  
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C--    Update Sea-Ice state :  
           IF ( compact.GT.0. _d 0 .AND. icFrac.EQ.0. _d 0) THEN  
              Tsrf(i,j,bi,bj)  = TFrzOce  
              Tice1(i,j,bi,bj) = TFrzOce  
              Tice2(i,j,bi,bj) = TFrzOce  
              Qice1(i,j,bi,bj) = qicen(1)  
              Qice2(i,j,bi,bj) = qicen(2)  
           ENDIF  
           iceheight(i,j,bi,bj) = hIce  
           snowheight(i,j,bi,bj)= hSnow  
328  C--    Net fluxes :  C--    Net fluxes :
329            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc  #ifdef ALLOW_ATM2D
330            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw            pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
331            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt            pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
332              sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
333            IF (dBug) WRITE(6,1010)'ThSI_FWD:-4- compact,hIc,hSn,Qnet=',  #else
334       &                        compact,hIce,hSnow,Qnet(i,j,bi,bj)            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
335  C--   - if esurp > 0 : end            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)
336              saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
337    #endif
338    
339    #ifdef ALLOW_DBUG_THSICE
340              IF (dBug(i,j,bi,bj)) WRITE(6,1010)
341         &         'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
342         &           iceFrac(i,j), iceHeight(i,j,bi,bj),
343         &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
344    #endif
345          ENDIF          ENDIF
346    
347          IF ( compact .GT. 0. _d 0 ) THEN          IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
348            iceMask(i,j,bi,bj)=compact       &    isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
349            IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0       &                  .AND.   iceFrac(i,j) .LE.0. _d 0
350            IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
351              iceMask(i,j,bi,bj)=iceFrac(i,j)
352              IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
353         &     snowAge(i,j,bi,bj) = 0. _d 0
354          ELSE          ELSE
355            iceMask(i,j,bi,bj)  = 0. _d 0            iceMask(i,j,bi,bj)  = 0. _d 0
356            iceHeight(i,j,bi,bj)= 0. _d 0            iceHeight(i,j,bi,bj)= 0. _d 0
357            snowHeight(i,j,bi,bj)=0. _d 0            snowHeight(i,j,bi,bj)=0. _d 0
358            snowAge(i,j,bi,bj)  = 0. _d 0            snowAge(i,j,bi,bj)  = 0. _d 0
359            Tsrf(i,j,bi,bj)     = oceTs            Tsrf(i,j,bi,bj)     = tOceMxL(i,j,bi,bj)
360            Tice1(i,j,bi,bj)    = 0. _d 0            Tice1(i,j,bi,bj)    = 0. _d 0
361            Tice2(i,j,bi,bj)    = 0. _d 0            Tice2(i,j,bi,bj)    = 0. _d 0
362            Qice1(i,j,bi,bj)    = 0. _d 0            Qice1(i,j,bi,bj)    = Lfresh
363            Qice2(i,j,bi,bj)    = 0. _d 0            Qice2(i,j,bi,bj)    = Lfresh
364          ENDIF          ENDIF
365           ENDDO
366          ENDDO
367    
368  C--     Return atmospheric fluxes in evpAtm & flxSW (same sign and units):  # ifdef ALLOW_AUTODIFF_TAMC
369          evpAtm(i,j) = frwAtm  CADJ STORE snowHeight(:,:,bi,bj) =
370          flxSW (i,j) = flxAtm(i,j)  CADJ &     comlev1_bibj,key=ticekey,byte=isbyte
371    # endif
372          DO j = jMin, jMax
373           DO i = iMin, iMax
374    C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
375            sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
376         &                         + iceHeight(i,j,bi,bj)*rhoi
377         &                        )*iceMask(i,j,bi,bj)
378         ENDDO         ENDDO
379        ENDDO        ENDDO
380    
381          IF ( thSIceAdvScheme.GT.0 ) THEN
382    C--   note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
383           DO j = jMin, jMax
384            DO i = iMin, iMax
385             IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
386              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
387              EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)
388              saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
389             ENDIF
390            ENDDO
391           ENDDO
392          ENDIF
393    
394    #ifdef ALLOW_BULK_FORCE
395          IF ( useBulkForce ) THEN
396            CALL BULKF_FLUX_ADJUST(
397         I                          bi, bj, iMin, iMax, jMin, jMax,
398         I                          isIceFree, myTime, myIter, myThid )
399          ENDIF
400    #endif /* ALLOW_BULK_FORCE */
401    
402  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
403  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
404    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22