/[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.16 by heimbach, Sun Apr 9 17:35:30 2006 UTC revision 1.21 by heimbach, Tue Apr 17 23:42:33 2007 UTC
# Line 9  C     !INTERFACE: Line 9  C     !INTERFACE:
9        SUBROUTINE THSICE_STEP_FWD(        SUBROUTINE THSICE_STEP_FWD(
10       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
11       I             prcAtm,       I             prcAtm,
      U             evpAtm, flxSW,  
12       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 30  C     === Global variables === Line 29  C     === Global variables ===
29  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
30  #include "THSICE_VARS.h"  #include "THSICE_VARS.h"
31  #include "THSICE_TAVE.h"  #include "THSICE_TAVE.h"
32    #include "THSICE_2DYN.h"
33    #ifdef ALLOW_AUTODIFF_TAMC
34    # include "tamc.h"
35    # include "tamc_keys.h"
36    #endif
37    
38          INTEGER siLo, siHi, sjLo, sjHi
39          PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
40          PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
41    
42  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
43  C     === Routine arguments ===  C     === Routine arguments ===
44    C- input:
45  C     bi,bj   :: tile indices  C     bi,bj   :: tile indices
46  C   iMin,iMax :: computation domain: 1rst index range  C   iMin,iMax :: computation domain: 1rst index range
47  C   jMin,jMax :: computation domain: 2nd  index range  C   jMin,jMax :: computation domain: 2nd  index range
 C- input:  
48  C     prcAtm  :: total precip from the atmosphere [kg/m2/s]  C     prcAtm  :: total precip from the atmosphere [kg/m2/s]
49  C     evpAtm  :: (Inp) evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)  C     myTime  :: current Time of simulation [s]
50  C     flxSW   :: (Inp) short-wave heat flux (+=down): downward comp. only  C     myIter  :: current Iteration number in simulation
51  C                      (part.1), becomes net SW flux into ocean (part.2).  C     myThid  :: my Thread Id number
52    C-- Use fluxes hold in commom blocks
53    C- input:
54    C     icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
55    C     icFlxAtm  :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
56    C     icFrwAtm  :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
57  C- output  C- output
58  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)
59  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.  
60        INTEGER bi,bj        INTEGER bi,bj
61        INTEGER iMin, iMax        INTEGER iMin, iMax
62        INTEGER jMin, jMax        INTEGER jMin, jMax
63        _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)  
64        _RL  myTime        _RL  myTime
65        INTEGER myIter        INTEGER myIter
66        INTEGER myThid        INTEGER myThid
# Line 61  CEOP Line 69  CEOP
69  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
70  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
71  C     === Local variables ===  C     === Local variables ===
72  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)  
73  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]
74  C     frw2oc    :: fresh-water flux from the ice to the ocean  C     frw2oc    :: fresh-water flux from the ice to the ocean
75  C     fsalt     :: mass salt flux to the ocean  C     fsalt     :: mass salt flux to the ocean
76  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]  C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
77  C     TFrzOce   :: sea-water freezing temperature [oC] (function of S)  C     tFrzOce   :: sea-water freezing temperature [oC] (function of S)
78  C     isIceFree :: true for ice-free grid-cell that remains ice-free  C     isIceFree :: true for ice-free grid-cell that remains ice-free
79        INTEGER i,j  C     ageFac    :: snow aging factor [1]
80        _RL snowPr  C     snowFac   :: snowing refreshing-age factor [units of 1/snowPr]
       _RL agingTime, ageFac  
       _RL albedo  
       _RL frwAtm  
       _RL flx2oc  
       _RL frw2oc  
       _RL fsalt  
       _RL TFrzOce, cphm, frzmltMxL  
       _RL Fbot, esurp  
       _RL opFrac, icFrac  
       _RL oceV2s, oceTs  
       _RL Tsf, Tice(nlyr), qicen(nlyr)  
       _RL tmpflx(0:2), tmpdTs  
       _RL hSnow  
 c  
       _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL compact(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL hIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
 c  
81        LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL     iceFrac  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL     flx2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84          _RL     frw2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85          _RL     fsalt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86          _RL     tFrzOce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RL     frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88          _RL ageFac
89          _RL snowFac
90          _RL cphm
91          _RL opFrac, icFrac
92  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
93        _RL tmpFac        _RL tmpFac
94  #endif  #endif
95          INTEGER i,j
96          LOGICAL dBugFlag
97    
98        LOGICAL dBug  C-    define grid-point location where to print debugging values
99    #include "THSICE_DEBUG.h"
100    
101   1010 FORMAT(A,1P4E14.6)   1010 FORMAT(A,1P4E14.6)
102        dBug = .FALSE.  
103  C-    Initialise flxAtm  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
104         DO j = 1-Oly, sNy+Oly  
105          DO i = 1-Olx, sNx+Olx  #ifdef ALLOW_AUTODIFF_TAMC
106            flxAtm(i,j) = 0.        act1 = bi - myBxLo(myThid)
107          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
108          act2 = bj - myByLo(myThid)
109          max2 = myByHi(myThid) - myByLo(myThid) + 1
110          act3 = myThid - 1
111          max3 = nTx*nTy
112          act4 = ikey_dynamics - 1
113          iicekey = (act1 + 1) + act2*max1
114         &                     + act3*max1*max2
115         &                     + act4*max1*max2*max3
116    #endif /* ALLOW_AUTODIFF_TAMC */
117    
118    C-    Initialise
119          dBugFlag = debugLevel.GE.debLevB
120          DO j = 1-OLy, sNy+OLy
121            DO i = 1-OLx, sNx+OLx
122            isIceFree(i,j) = .FALSE.            isIceFree(i,j) = .FALSE.
123              saltFlux(i,j,bi,bj) = 0. _d 0
124  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
125            compact(i,j) = 0.            iceFrac(i,j) = 0.
           hIce(i,j) = 0.  
126  #endif  #endif
127          ENDDO          ENDDO
128         ENDDO        ENDDO
   
       IF ( fluidIsWater ) THEN  
        DO j = jMin, jMax  
         DO i = iMin, iMax  
 c        dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )  
129    
130  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|        ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
131  C    part.1 : ice-covered fraction ;        snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
 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(i,j)    = 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)  
           IF ( dBug ) THEN  
            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj  
            WRITE(6,1010) 'ThSI_FWD:-0- iceMask, hIc, hSn, Tsf  =',  
      &                                 icFrac, hIce,hSnow,Tsf  
           ENDIF  
   
           CALL THSICE_ALBEDO(  
      I               hIce(i,j), hSnow, Tsf, snowAge(i,j,bi,bj),  
      O               albedo,  
      I               myThid )  
           flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo)  
           siceAlb(i,j,bi,bj) = albedo  
   
           CALL THSICE_SOLVE4TEMP(  
      I          useBulkForce, tmpflx, TFrzOce, hIce(i,j), 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)  
132    
133  #ifdef SHORTWAVE_HEATING  #ifdef ALLOW_AUTODIFF_TAMC
134  C--    Update Fluxes :  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
           opFrac= 1. _d 0-icFrac  
           Qsw(i,j,bi,bj)=-icFrac*flxSW(i,j) +opFrac*Qsw(i,j,bi,bj)  
135  #endif  #endif
136  C--    Update Sea-Ice state :        DO j = jMin, jMax
137            Tsrf(i,j,bi,bj) =Tsf         DO i = iMin, iMax
138            Tice1(i,j,bi,bj)=Tice(1)          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
139            Tice2(i,j,bi,bj)=Tice(2)  C--  Snow aging :
140            Qice1(i,j,bi,bj)=qicen(1)            snowAge(i,j,bi,bj) = thSIce_deltaT
141            Qice2(i,j,bi,bj)=qicen(2)       &                       + snowAge(i,j,bi,bj)*ageFac
142            IF ( dBug ) THEN            IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
143             WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), TFrzOce =',       &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
144       &                              Tsf, Tice, TFrzOce       &          * EXP( - snowFac*snowPrc(i,j,bi,bj) )
145             WRITE(6,1010) 'ThSI_FWD: sHeat,fxCndBt, fxAtm,evAtm=',  c    &          * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
146       &                  sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),  c    &                  /hNewSnowAge )
147       &                  flxAtm(i,j), evpAtm(i,j)  C-------
148            ENDIF  C note: Any flux of mass (here fresh water) that enter or leave the system
149           ENDIF  C       with a non zero energy HAS TO be counted: add snow precip.
150          ENDDO            icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
151         &                        - Lfresh*snowPrc(i,j,bi,bj)
152    C--
153            ENDIF
154         ENDDO         ENDDO
155        ENDIF        ENDDO
       dBug = .FALSE.  
156    
157  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
158        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
# Line 198  C    part.2 : ice-covered fraction ; Line 176  C    part.2 : ice-covered fraction ;
176  C     change in ice/snow thickness and ice-fraction  C     change in ice/snow thickness and ice-fraction
177  C     note: can only reduce the ice-fraction but not increase it.  C     note: can only reduce the ice-fraction but not increase it.
178  C-------  C-------
       agingTime = 50. _d 0 * 86400. _d 0  
       ageFac = 1. _d 0 - thSIce_deltaT/agingTime  
179        DO j = jMin, jMax        DO j = jMin, jMax
180         DO i = iMin, iMax         DO i = iMin, iMax
 c       dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )  
181    
182          TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)          tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
         oceTs   = tOceMxL(i,j,bi,bj)  
183          cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)          cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)
184          frzmltMxL = (TFrzOce-oceTs)*cphm/ocean_deltaT          frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
185         &                 * cphm/ocean_deltaT
186          Fbot   = 0. _d 0          iceFrac(i,j) = iceMask(i,j,bi,bj)
187          saltFlux(i,j,bi,bj) = 0. _d 0          flx2oc(i,j)  = icFlxSW(i,j,bi,bj)
         compact(i,j)= iceMask(i,j,bi,bj)  
188  C-------  C-------
189          IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact(i,j).GT.0.) ) THEN  #ifdef ALLOW_DBUG_THSICE
190            IF ( dBug(i,j,bi,bj) ) THEN
191             IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
192            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj            WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
193            WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',            WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
194       &                  compact(i,j), iceHeight(i,j,bi,bj),       &                  iceFrac(i,j), iceHeight(i,j,bi,bj),
195       &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)       &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
196            WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL,Qnet=',            WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
197       &                     oceTs, TFrzOce, frzmltMxL,Qnet(i,j,bi,bj)       &                     tOceMxL(i,j,bi,bj), tFrzOce(i,j),
198         &                     frzmltMxL(i,j), Qnet(i,j,bi,bj)
199             ENDIF
200             IF (iceFrac(i,j).GT.0.)
201         &    WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
202         &      iceFrac(i,j), icFlxAtm(i,j,bi,bj),
203         &      icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
204          ENDIF          ENDIF
205  C-------  #endif
206          IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN         ENDDO
207          ENDDO
           oceV2s  = v2ocMxL(i,j,bi,bj)  
           snowPr  = snowPrc(i,j,bi,bj)  
           hIce(i,j)    = 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(i,j), hIce(i,j), 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)  
208    
209  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  #ifdef ALLOW_AUTODIFF_TAMC
210        IF (dBug) WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',  CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
211       &  iceMask(i,j,bi,bj),flxAtm(i,j),evpAtm(i,j),-Lfresh*snowPr  #endif
       IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc=',  
      &   compact(i,j),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(i,j), hIce(i,j), 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-|--+----|  
212    
213  C--    Update Sea-Ice state :        CALL THSICE_CALC_THICKN(
214  c         iceMask(i,j,bi,bj)=compact       I          bi, bj, siLo, siHi, sjLo, sjHi,
215            iceHeight(i,j,bi,bj) = hIce(i,j)       I          iMin,iMax, jMin,jMax, dBugFlag,
216            snowHeight(i,j,bi,bj)= hSnow       I          iceMask(siLo,sjLo,bi,bj), tFrzOce,
217            Tsrf(i,j,bi,bj) =Tsf       I          tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
218            Qice1(i,j,bi,bj)=qicen(1)       I          snowPrc(siLo,sjLo,bi,bj), prcAtm,
219            Qice2(i,j,bi,bj)=qicen(2)       I          sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
220         U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
221         U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
222         U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
223         U          icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
224         O          frw2oc, fsalt,
225         I          myTime, myIter, myThid )
226    
227  C--    Net fluxes :  C--    Net fluxes :
228            frw2oc = frw2oc + (prcAtm(i,j)-snowPrc(i,j,bi,bj))        DO j = jMin, jMax
229           DO i = iMin, iMax
230    #ifdef ALLOW_AUTODIFF_TAMC
231              ikey_1 = i
232         &         + sNx*(j-1)
233         &         + sNx*sNy*act1
234         &         + sNx*sNy*max1*act2
235         &         + sNx*sNy*max1*max2*act3
236         &         + sNx*sNy*max1*max2*max3*act4
237    #endif /* ALLOW_AUTODIFF_TAMC */
238    C--
239    #ifdef ALLOW_AUTODIFF_TAMC
240    CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
241    #endif
242            IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
243  C-     weighted average net fluxes:  C-     weighted average net fluxes:
244    #ifdef ALLOW_AUTODIFF_TAMC
245    CADJ STORE  fsalt(i,j) = comlev1_thsice_1, key=ikey_1
246    CADJ STORE  flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
247    CADJ STORE  frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
248    CADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
249    #endif
250            icFrac = iceMask(i,j,bi,bj)            icFrac = iceMask(i,j,bi,bj)
251            opFrac= 1. _d 0-icFrac            opFrac= 1. _d 0-icFrac
252            flxAtm(i,j) = icFrac*flxAtm(i,j) - opFrac*Qnet(i,j,bi,bj)            icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
253            frwAtm =     icFrac*frwAtm + opFrac*rhofw*EmPmR(i,j,bi,bj)       &                        - opFrac*Qnet(i,j,bi,bj)
254            Qnet(i,j,bi,bj)=-icFrac*flx2oc +opFrac*Qnet(i,j,bi,bj)            icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
255            EmPmR(i,j,bi,bj)=-icFrac*frw2oc/rhofw+opFrac*EmPmR(i,j,bi,bj)       &                        + opFrac*rhofw*EmPmR(i,j,bi,bj)
256            saltFlux(i,j,bi,bj)=-icFrac*fsalt            Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
257              EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
258            IF (dBug) WRITE(6,1010)       &                    +  opFrac*EmPmR(i,j,bi,bj)
259       &          'ThSI_FWD:-3- compact, hIc, hSn, Qnet =',            saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
260       &                        compact(i,j),hIce(i,j),hSnow,  
261       &                        Qnet(i,j,bi,bj)  #ifdef ALLOW_DBUG_THSICE
262              IF (dBug(i,j,bi,bj)) WRITE(6,1010)
263         &          'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
264         &           iceFrac(i,j), iceHeight(i,j,bi,bj),
265         &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
266    #endif
267    
268          ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN          ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
269            flxAtm(i,j) =  -Qnet(i,j,bi,bj)            icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
270            frwAtm = rhofw*EmPmR(i,j,bi,bj)            icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
271          ELSE          ELSE
272            flxAtm(i,j) = 0. _d 0            icFlxAtm(i,j,bi,bj) = 0. _d 0
273            frwAtm      = 0. _d 0            icFrwAtm(i,j,bi,bj) = 0. _d 0
274          ENDIF          ENDIF
275           ENDDO
276          ENDDO
277    
278  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
279  C    part.3 : freezing of sea-water  C    part.3 : freezing of sea-water
280  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
281  C-------  C-------
282  c       compact(i,j) = iceMask(i,j,bi,bj)        CALL THSICE_EXTEND(
283          hIce(i,j)   = iceHeight(i,j,bi,bj)       I          bi, bj, siLo, siHi, sjLo, sjHi,
284          hSnow  = snowHeight(i,j,bi,bj)       I          iMin,iMax, jMin,jMax, dBugFlag,
285         I          frzmltMxL, tFrzOce,
286          esurp  = frzmltMxL - Fbot*iceMask(i,j,bi,bj)       I          tOceMxL(siLo,sjLo,bi,bj),
287          IF (esurp.GT.0. _d 0) THEN       U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
288            icFrac = compact(i,j)       U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
289            qicen(1)= Qice1(i,j,bi,bj)       U          Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
290            qicen(2)= Qice2(i,j,bi,bj)       U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
291            CALL THSICE_EXTEND(       O          flx2oc, frw2oc, fsalt,
292       I               esurp, TFrzOce,       I          myTime, myIter, myThid )
293       U               oceTs, compact(i,j), hIce(i,j), hSnow, qicen,  
294       O               flx2oc, frw2oc, fsalt,  #ifdef ALLOW_AUTODIFF_TAMC
295       I               dBug, myThid )  CADJ STORE snowHeight(:,:,bi,bj) =
296  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
297        IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc='  #endif
298       &                 ,compact(i,j),flx2oc,fsalt,frw2oc        DO j = jMin, jMax
299  #ifdef CHECK_ENERGY_CONSERV         DO i = iMin, iMax
300            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(i,j), hIce(i,j), 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(i,j).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(i,j)  
           snowHeight(i,j,bi,bj)= hSnow  
301  C--    Net fluxes :  C--    Net fluxes :
302            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc            Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
303            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw            EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
304            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt            saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
305    
306            IF (dBug) WRITE(6,1010)  #ifdef ALLOW_DBUG_THSICE
307       &         'ThSI_FWD:-4- compact, hIc, hSn, Qnet =',            IF (dBug(i,j,bi,bj)) WRITE(6,1010)
308       &         compact(i,j),hIce(i,j),hSnow,       &         'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
309       &         Qnet(i,j,bi,bj)       &           iceFrac(i,j), iceHeight(i,j,bi,bj),
310  C--   - if esurp > 0 : end       &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
311    #endif
312          ENDIF          ENDIF
313    
314          IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )          IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
315       &    isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0       &    isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
316       &                       .AND.   compact(i,j) .LE.0. _d 0       &                  .AND.   iceFrac(i,j) .LE.0. _d 0
317          IF ( compact(i,j) .GT. 0. _d 0 ) THEN          IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
318            iceMask(i,j,bi,bj)=compact(i,j)            iceMask(i,j,bi,bj)=iceFrac(i,j)
319            IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0            IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
320         &     snowAge(i,j,bi,bj) = 0. _d 0
321          ELSE          ELSE
322            iceMask(i,j,bi,bj)  = 0. _d 0            iceMask(i,j,bi,bj)  = 0. _d 0
323            iceHeight(i,j,bi,bj)= 0. _d 0            iceHeight(i,j,bi,bj)= 0. _d 0
324            snowHeight(i,j,bi,bj)=0. _d 0            snowHeight(i,j,bi,bj)=0. _d 0
325            snowAge(i,j,bi,bj)  = 0. _d 0            snowAge(i,j,bi,bj)  = 0. _d 0
326            Tsrf(i,j,bi,bj)     = oceTs            Tsrf(i,j,bi,bj)     = tOceMxL(i,j,bi,bj)
327            Tice1(i,j,bi,bj)    = 0. _d 0            Tice1(i,j,bi,bj)    = 0. _d 0
328            Tice2(i,j,bi,bj)    = 0. _d 0            Tice2(i,j,bi,bj)    = 0. _d 0
329            Qice1(i,j,bi,bj)    = 0. _d 0            Qice1(i,j,bi,bj)    = Lfresh
330            Qice2(i,j,bi,bj)    = 0. _d 0            Qice2(i,j,bi,bj)    = Lfresh
331          ENDIF          ENDIF
332           ENDDO
333  C--     Return atmospheric fluxes in evpAtm & flxSW (same sign and units):        ENDDO
         evpAtm(i,j) = frwAtm  
         flxSW (i,j) = flxAtm(i,j)  
334    
335  #ifdef ATMOSPHERIC_LOADING  #ifdef ATMOSPHERIC_LOADING
336    # ifdef ALLOW_AUTODIFF_TAMC
337    CADJ STORE snowHeight(:,:,bi,bj) =
338    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
339    # endif
340          DO j = jMin, jMax
341           DO i = iMin, iMax
342  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)  C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
343          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos          sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
344       &                         + iceHeight(i,j,bi,bj)*rhoi       &                         + iceHeight(i,j,bi,bj)*rhoi
345       &                        )*iceMask(i,j,bi,bj)       &                        )*iceMask(i,j,bi,bj)
346           ENDDO
347          ENDDO
348  #endif  #endif
349    
350          IF ( thSIceAdvScheme.GT.0 ) THEN
351    C--   note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
352           DO j = jMin, jMax
353            DO i = iMin, iMax
354             IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
355              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
356              EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)/rhofw
357              saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
358             ENDIF
359            ENDDO
360         ENDDO         ENDDO
361        ENDDO        ENDIF
362    
363  #ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_BULK_FORCE
364        IF ( useBulkForce ) THEN        IF ( useBulkForce ) THEN
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22