/[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.15 by heimbach, Tue Apr 17 23:42:33 2007 UTC revision 1.16 by jmc, Sun Apr 29 23:48:44 2007 UTC
# Line 212  C---+----1----+----2----+----3----+----4 Line 212  C---+----1----+----2----+----3----+----4
212            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
213       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
214  #endif  #endif
215            IF ( hi.LT.himin ) THEN            IF ( hi.LT.hIceMin ) THEN
216  C If hi < himin, melt the ice.  C If hi < hIceMin, melt the ice.
217               STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'               STOP 'THSICE_SOLVE4TEMP: should not enter if hi<hIceMin'
218            ENDIF            ENDIF
219    
220  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 224  C fractional snow cover Line 224  C fractional snow cover
224        IF (hs .GT. 0. _d 0) frsnow = 1. _d 0        IF (hs .GT. 0. _d 0) frsnow = 1. _d 0
225    
226  C Compute SW flux absorbed at surface and penetrating to layer 1.  C Compute SW flux absorbed at surface and penetrating to layer 1.
227        fswpen  = netSW * (1. _d 0 - frsnow) * i0        fswpen  = netSW * (1. _d 0 - frsnow) * i0swFrac
228        fswocn = fswpen * exp(-ksolar*hi)        fswocn = fswpen * exp(-ksolar*hi)
229        fswint = fswpen - fswocn        fswint = fswpen - fswocn
230    
# Line 232  C Compute SW flux absorbed at surface an Line 232  C Compute SW flux absorbed at surface an
232    
233  C Compute conductivity terms at layer interfaces.  C Compute conductivity terms at layer interfaces.
234    
235        k12 = 4. _d 0*kice*ksnow / (ksnow*hi + 4. _d 0*kice*hs)        k12 = 4. _d 0*kIce*kSnow / (kSnow*hi + 4. _d 0*kIce*hs)
236        k32 = 2. _d 0*kice  / hi        k32 = 2. _d 0*kIce  / hi
237    
238  C compute ice temperatures  C compute ice temperatures
239        a1 = cpice        a1 = cpIce
240        b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh        b1 = qicen(1) + (cpWater-cpIce )*Tmlt1 - Lfresh
241        c1 = Lfresh * Tmlt1        c1 = Lfresh * Tmlt1
242        Tice(1) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1        Tice(1) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
243        Tice(2) = (Lfresh-qicen(2)) / cpice        Tice(2) = (Lfresh-qicen(2)) / cpIce
244    
245        IF (Tice(1).GT.0. _d 0 ) THEN        IF (Tice(1).GT.0. _d 0 ) THEN
246         WRITE (standardMessageUnit,'(A,I12,1PE14.6)')         WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
# Line 261  C compute ice temperatures Line 261  C compute ice temperatures
261    
262  C Compute coefficients used in quadratic formula.  C Compute coefficients used in quadratic formula.
263    
264        a10 = rhoi*cpice *hi/(2. _d 0*dt) +        a10 = rhoi*cpIce *hi/(2. _d 0*dt) +
265       &      k32 * (4. _d 0*dt*k32 + rhoi*cpice *hi)       &      k32 * (4. _d 0*dt*k32 + rhoi*cpIce *hi)
266       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi)       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi)
267        b10 = -hi*        b10 = -hi*
268       &      (rhoi*cpice*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))       &      (rhoi*cpIce*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))
269       &       /(2. _d 0*dt)       &       /(2. _d 0*dt)
270       &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpice *hi*Tice(2))       &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpIce *hi*Tice(2))
271       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi) - fswint       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi) - fswint
272        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)
273    
274  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 410  CADJ STORE  Tice(:)      = comlev1_thsic Line 410  CADJ STORE  Tice(:)      = comlev1_thsic
410  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1
411  #endif  #endif
412        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)
413       &        + rhoi*cpice *hi*Tice(2))       &        + rhoi*cpIce *hi*Tice(2))
414       &         /(6. _d 0*dt*k32 + rhoi*cpice *hi)       &         /(6. _d 0*dt*k32 + rhoi*cpIce *hi)
415  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
416        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
417       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice
# Line 420  CADJ STORE  df0dt        = comlev1_thsic Line 420  CADJ STORE  df0dt        = comlev1_thsic
420  C Compute final flux values at surfaces.  C Compute final flux values at surfaces.
421    
422        fct    = k12*(Tsf-Tice(1))        fct    = k12*(Tsf-Tice(1))
423        flxCnB(i,j) = 4. _d 0*kice *(Tice(2)-Tf)/hi        flxCnB(i,j) = 4. _d 0*kIce *(Tice(2)-Tf)/hi
424        flx0   = flx0 + df0dT*dTsf        flx0   = flx0 + df0dT*dTsf
425        IF ( useBlkFlx ) THEN        IF ( useBlkFlx ) THEN
426  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
# Line 447  C-    SW flux at sea-ice base left to th Line 447  C-    SW flux at sea-ice base left to th
447    
448  C Compute new enthalpy for each layer.  C Compute new enthalpy for each layer.
449    
450        qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1))        qicen(1) = -cpWater*Tmlt1 + cpIce *(Tmlt1-Tice(1))
451       &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))       &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))
452        qicen(2) = -cpice *Tice(2) + Lfresh        qicen(2) = -cpIce *Tice(2) + Lfresh
453    
454  C Make sure internal ice temperatures do not exceed Tmlt.  C Make sure internal ice temperatures do not exceed Tmlt.
455  C (This should not happen for reasonable values of i0.)  C (This should not happen for reasonable values of i0swFrac)
456    
457        IF (Tice(1) .GE. Tmlt1) THEN        IF (Tice(1) .GE. Tmlt1) THEN
458          WRITE (6,'(A,2I4,2I3,1P2E14.6)')          WRITE (6,'(A,2I4,2I3,1P2E14.6)')
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22