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

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

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

revision 1.19 by jmc, Tue Mar 16 00:23:59 2010 UTC revision 1.20 by jmc, Thu Sep 30 16:05:19 2010 UTC
# Line 140  C     fswpen       :: SW penetrating ben Line 140  C     fswpen       :: SW penetrating ben
140  C     fswdn        :: SW absorbed at surface (W m-2)  C     fswdn        :: SW absorbed at surface (W m-2)
141  C     fswint       :: SW absorbed in ice (W m-2)  C     fswint       :: SW absorbed in ice (W m-2)
142  C     fswocn       :: SW passed through ice to ocean (W m-2)  C     fswocn       :: SW passed through ice to ocean (W m-2)
143  C     flxExceptSw  :: net surface heat flux, except short-wave (W/m2)  C     flx0exSW     :: net surface heat flux over melting snow/ice, except short-wave.
144  C     evap         :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)  C     flxTexSW     :: net surface heat flux, except short-wave (W/m2)
145    C     evap0        :: evaporation over melting snow/ice [kg/m2/s] (>0 if evaporate)
146    C     evapT        :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
147  C     dEvdT        :: derivative of evap. with respect to Tsf [kg/m2/s/K]  C     dEvdT        :: derivative of evap. with respect to Tsf [kg/m2/s/K]
148  C     flx0         :: net surf heat flux, from Atmos. to sea-ice (W m-2)  C     flxNet       :: net surf heat flux (+=down), from Atmos. to sea-ice (W m-2)
149  C     fct          :: heat conducted to top surface  C     fct          :: heat conducted to top surface
150  C     df0dT        :: deriv of flx0 wrt Tsf (W m-2 deg-1)  C     dFlxdT       :: deriv of flxNet wrt Tsf (W m-2 deg-1)
151  C     k12, k32     :: thermal conductivity terms  C     k12, k32     :: thermal conductivity terms
152  C     a10, b10     :: coefficients in quadratic eqn for T1  C     a10, b10     :: coefficients in quadratic eqn for T1
153  C     a1, b1, c1   :: coefficients in quadratic eqn for T1  C     a1, b1, c1   :: coefficients in quadratic eqn for T1
# Line 158  C     dt           :: timestep Line 160  C     dt           :: timestep
160        _RL  fswdn        _RL  fswdn
161        _RL  fswint        _RL  fswint
162        _RL  fswocn        _RL  fswocn
163        _RL  flxExceptSw        _RL  flx0exSW
164        _RL  evap, dEvdT        _RL  flxTexSW
165        _RL  flx0        _RL  evap0, evapT, dEvdT
166          _RL  flxNet
167        _RL  fct        _RL  fct
168        _RL  df0dT        _RL  dFlxdT
169        _RL  k12, k32        _RL  k12, k32
170        _RL  a10, b10        _RL  a10, b10
171        _RL  a1, b1, c1        _RL  a1, b1, c1
172  c     _RL  Tsf_start  c     _RL  Tsf_start
173        _RL  dt        _RL  dt
174        _RL  recip_dhSnowLin        _RL  recip_dhSnowLin
       INTEGER iceornot  
175        LOGICAL useBlkFlx        LOGICAL useBlkFlx
176    
177  C-    define grid-point location where to print debugging values  C-    define grid-point location where to print debugging values
# Line 211  C---+----1----+----2----+----3----+----4 Line 213  C---+----1----+----2----+----3----+----4
213  C--  C--
214  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
215  CADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1  CADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1
216  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1  CADJ STORE  dFlxdT        = comlev1_thsice_1, key=ikey_1
217  CADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1  CADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1
218  CADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1  CADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1
219  CADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1  CADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1
# Line 223  CADJ STORE  tsrf(i,j)    = comlev1_thsic Line 225  CADJ STORE  tsrf(i,j)    = comlev1_thsic
225            hs      = hSnow(i,j)            hs      = hSnow(i,j)
226            Tf      = tFrz(i,j)            Tf      = tFrz(i,j)
227            netSW   = flxSW(i,j)            netSW   = flxSW(i,j)
           Tsf     = tSrf(i,j)  
228            qicen(1)= qIc1(i,j)            qicen(1)= qIc1(i,j)
229            qicen(2)= qIc2(i,j)            qicen(2)= qIc2(i,j)
230  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 282  C compute ice temperatures Line 283  C compute ice temperatures
283        ENDIF        ENDIF
284  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
285        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
286       &  'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice       &  'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf(i,j),Tice
287  #endif  #endif
288    
289  C Compute coefficients used in quadratic formula.  C Compute coefficients used in quadratic formula.
# Line 298  C Compute coefficients used in quadratic Line 299  C Compute coefficients used in quadratic
299        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)
300    
301  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
302    C get surface fluxes over melting surface
303              IF ( useBlkFlx ) THEN
304               Tsf = 0.
305               IF ( useEXF ) THEN
306                CALL THSICE_GET_EXF(
307         I                            hSnow(i,j), Tsf,
308         O                            flx0exSW, dFlxdT, evap0, dEvdT,
309         I                            i,j,bi,bj,myThid )
310    C could add this "ifdef" to hide THSICE_GET_BULKF from TAF
311    c#ifdef ALLOW_BULK_FORCE
312               ELSEIF ( useBulkForce ) THEN
313                CALL THSICE_GET_BULKF(
314         I                            hSnow(i,j), Tsf,
315         O                            flx0exSW, dFlxdT, evap0, dEvdT,
316         I                            i,j,bi,bj,myThid )
317    c#endif /* ALLOW_BULK_FORCE */
318               ENDIF
319              ELSE
320               flx0exSW = flxExSW(i,j,0)
321              ENDIF
322    
323    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
324  C Compute new surface and internal temperatures; iterate until  C Compute new surface and internal temperatures; iterate until
325  C Tsfc converges.  C Tsfc converges.
326    
# Line 306  C Tsfc converges. Line 329  C Tsfc converges.
329        ELSE        ELSE
330          iterMax = 1          iterMax = 1
331        ENDIF        ENDIF
332          Tsf  = tSrf(i,j)
333        dTsf = Terrmax        dTsf = Terrmax
334    
335  C ----- begin iteration  -----  C ----- begin iteration  -----
# Line 318  C ----- begin iteration  ----- Line 342  C ----- begin iteration  -----
342  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
343  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3
344  CADJ STORE dtsf         = comlev1_thsice_3, key=ikey_3  CADJ STORE dtsf         = comlev1_thsice_3, key=ikey_3
345  CADJ STORE df0dt        = comlev1_thsice_3, key=ikey_3  CADJ STORE dFlxdT        = comlev1_thsice_3, key=ikey_3
346  CADJ STORE flxexceptsw  = comlev1_thsice_3, key=ikey_3  CADJ STORE flxexceptsw  = comlev1_thsice_3, key=ikey_3
347  #endif  #endif
348         IF ( ABS(dTsf).GE.Terrmax ) THEN         IF ( ABS(dTsf).GE.Terrmax ) THEN
# Line 331  c        Tsf_start = Tsf Line 355  c        Tsf_start = Tsf
355  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3
356  #endif  #endif
357  C Compute top surface flux.  C Compute top surface flux.
358           IF (hs.GT.3. _d -1) THEN           IF ( useEXF ) THEN
               iceornot=2  
          ELSE  
               iceornot=1  
          ENDIF  
          IF ( useBulkForce ) THEN  
           CALL THSICE_GET_BULKF(  
      I                         iceornot, Tsf,  
      O                         flxExceptSw, df0dT, evap, dEvdT,  
      I                         i,j,bi,bj,myThid )  
          ELSEIF ( useEXF ) THEN  
359            CALL THSICE_GET_EXF  (            CALL THSICE_GET_EXF  (
360       I                         iceornot, Tsf,       I                         hs, Tsf,
361       O                         flxExceptSw, df0dT, evap, dEvdT,       O                         flxTexSW, dFlxdT, evapT, dEvdT,
362         I                         i,j,bi,bj,myThid )
363    C could add this "ifdef" to hide THSICE_GET_BULKF from TAF
364    c#ifdef ALLOW_BULK_FORCE
365             ELSEIF ( useBulkForce ) THEN
366              CALL THSICE_GET_BULKF(
367         I                         hs, Tsf,
368         O                         flxTexSW, dFlxdT, evapT, dEvdT,
369       I                         i,j,bi,bj,myThid )       I                         i,j,bi,bj,myThid )
370    c#endif /* ALLOW_BULK_FORCE */
371           ENDIF           ENDIF
372          ELSE          ELSE
373           flxExceptSw = flxExSW(i,j,1)           flxTexSW = flxExSW(i,j,1)
374           df0dT = flxExSW(i,j,2)           dFlxdT = flxExSW(i,j,2)
375          ENDIF          ENDIF
376          flx0 = fswdn + flxExceptSw          flxNet = fswdn + flxTexSW
377  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
378        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
379       &  'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT       &  'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',
380         &                 flxNet, dFlxdT, k12, k12-dFlxdT
381  #endif  #endif
382    
383  C Compute new top layer and surface temperatures.  C Compute new top layer and surface temperatures.
# Line 364  C with different coefficients. Line 387  C with different coefficients.
387  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
388  CADJ STORE tsf  = comlev1_thsice_3, key=ikey_3  CADJ STORE tsf  = comlev1_thsice_3, key=ikey_3
389  #endif  #endif
390           a1 = a10 - k12*df0dT / (k12-df0dT)           a1 = a10 - k12*dFlxdT / (k12-dFlxdT)
391           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)           b1 = b10 - k12*(flxNet-dFlxdT*Tsf) / (k12-dFlxdT)
392           Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)           Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
393           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)           dTsf = (flxNet + k12*(Tice(1)-Tsf)) / (k12-dFlxdT)
394           Tsf = Tsf + dTsf           Tsf = Tsf + dTsf
395           IF (Tsf .GT. 0. _d 0) THEN           IF (Tsf .GT. 0. _d 0) THEN
396  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
# Line 380  C      note: b1 = b10 - k12*Tf0 Line 403  C      note: b1 = b10 - k12*Tf0
403              Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)              Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
404              Tsf = 0. _d 0              Tsf = 0. _d 0
405             IF ( useBlkFlx ) THEN             IF ( useBlkFlx ) THEN
406              IF (hs.GT.3. _d -1) THEN              flxTexSW = flx0exSW
407                   iceornot=2              evapT = evap0
             ELSE  
                  iceornot=1  
             ENDIF  
             IF ( useBulkForce ) THEN  
              CALL THSICE_GET_BULKF(  
      I                            iceornot, Tsf,  
      O                            flxExceptSw, df0dT, evap, dEvdT,  
      I                            i,j,bi,bj,myThid )  
             ELSEIF ( useEXF ) THEN  
              CALL THSICE_GET_EXF  (  
      I                            iceornot, Tsf,  
      O                            flxExceptSw, df0dT, evap, dEvdT,  
      I                            i,j,bi,bj,myThid )  
             ENDIF  
408              dTsf = 0. _d 0              dTsf = 0. _d 0
409             ELSE             ELSE
410              flxExceptSw = flxExSW(i,j,0)              flxTexSW = flxExSW(i,j,0)
411              dTsf = 1000.              dTsf = 1000.
412              df0dT = 0.              dFlxdT = 0.
413             ENDIF             ENDIF
414             flx0 = fswdn + flxExceptSw             flxNet = fswdn + flxTexSW
415           ENDIF           ENDIF
416    
417  C Check for convergence.  If no convergence, then repeat.  C Check for convergence.  If no convergence, then repeat.
# Line 422  C If no convergence, then repeat. Line 431  C If no convergence, then repeat.
431       &                 ' BB: not converge: i,j,it,hi=',i,j,bi,bj,       &                 ' BB: not converge: i,j,it,hi=',i,j,bi,bj,
432       &                                                    myIter,hi       &                                                    myIter,hi
433              WRITE (6,*) 'BB: not converge: Tsf, dTsf=', Tsf,dTsf              WRITE (6,*) 'BB: not converge: Tsf, dTsf=', Tsf,dTsf
434              WRITE (6,*) 'BB: not converge: flx0,dfdT=',flx0,df0dT              WRITE (6,*) 'BB: not converge: flxNet,dFlxT=',flxNet,dFlxdT
435              IF (Tsf.LT.-70. _d 0) STOP              IF (Tsf.LT.-70. _d 0) STOP
436           ENDIF           ENDIF
437    
# Line 434  C Compute new bottom layer temperature. Line 443  C Compute new bottom layer temperature.
443    
444  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
445  CADJ STORE  Tice(:)      = comlev1_thsice_1, key=ikey_1  CADJ STORE  Tice(:)      = comlev1_thsice_1, key=ikey_1
446  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1  CADJ STORE  dFlxdT        = comlev1_thsice_1, key=ikey_1
447  #endif  #endif
448        Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)        Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)
449       &        + rhoi*cpIce *hi*Tice(2))       &        + rhoi*cpIce *hi*Tice(2))
# Line 448  C Compute final flux values at surfaces. Line 457  C Compute final flux values at surfaces.
457    
458        fct    = k12*(Tsf-Tice(1))        fct    = k12*(Tsf-Tice(1))
459        flxCnB(i,j) = 4. _d 0*kIce *(Tice(2)-Tf)/hi        flxCnB(i,j) = 4. _d 0*kIce *(Tice(2)-Tf)/hi
460        flx0   = flx0 + df0dT*dTsf        flxNet   = flxNet + dFlxdT*dTsf
461        IF ( useBlkFlx ) THEN        IF ( useBlkFlx ) THEN
462  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated
463          evpAtm(i,j) = evap + dEvdT*dTsf          evpAtm(i,j) = evapT + dEvdT*dTsf
464        ELSE        ELSE
465  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
466          evpAtm(i,j) = 0.          evpAtm(i,j) = 0.
467        ENDIF        ENDIF
468  C-    energy flux to Atmos: use net short-wave flux at surf. and  C-    energy flux to Atmos: use net short-wave flux at surf. and
469  C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)  C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)
470        flxAtm(i,j) = netSW + flxExceptSw        flxAtm(i,j) = netSW + flxTexSW
471       &                    + df0dT*dTsf + evpAtm(i,j)*Lfresh       &                    + dFlxdT*dTsf + evpAtm(i,j)*Lfresh
472  C-    excess of energy @ surface (used for surface melting):  C-    excess of energy @ surface (used for surface melting):
473        sHeat(i,j) = flx0 - fct        sHeat(i,j) = flxNet - fct
474    
475  C-    SW flux at sea-ice base left to the ocean  C-    SW flux at sea-ice base left to the ocean
476        flxSW(i,j) = fswocn        flxSW(i,j) = fswocn
477    
478  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
479        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
480       &  'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',       &  'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
481       &                 flx0,fct,flx0-fct,flxCnB(i,j)       &                 flxNet,fct,flxNet-fct,flxCnB(i,j)
482  #endif  #endif
483    
484  C Compute new enthalpy for each layer.  C Compute new enthalpy for each layer.
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22