/[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.21 by jmc, Sun Oct 10 20:15:06 2010 UTC revision 1.22 by mlosch, Fri Oct 15 15:15:31 2010 UTC
# Line 146  C     dt           :: timestep Line 146  C     dt           :: timestep
146        LOGICAL iterate4Tsf        LOGICAL iterate4Tsf
147        LOGICAL iceFlag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        LOGICAL iceFlag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        INTEGER i,j        INTEGER i,j
149        INTEGER k, iterMax        INTEGER k, iterMax, ii, jj, icount
150        _RL  frsnow        _RL  frsnow
151        _RL  fswpen        _RL  fswpen
152        _RL  fswdn        _RL  fswdn
# Line 200  c     act4 = ikey_dynamics - 1 Line 200  c     act4 = ikey_dynamics - 1
200        ENDIF        ENDIF
201    
202        iterate4Tsf = .FALSE.        iterate4Tsf = .FALSE.
203          icount = 0
204    C    
205        DO j = jMin, jMax        DO j = jMin, jMax
206         DO i = iMin, iMax         DO i = iMin, iMax
207  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 227  cCADJ STORE  tsrf(i,j)    = comlev1_thsi Line 229  cCADJ STORE  tsrf(i,j)    = comlev1_thsi
229       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
230  #endif  #endif
231            IF ( hIce(i,j).LT.hIceMin ) THEN            IF ( hIce(i,j).LT.hIceMin ) THEN
232  C         if hi < hIceMin, melt the ice.  C     if hi < hIceMin, melt the ice.
233               STOP 'THSICE_SOLVE4TEMP: should not enter if hi<hIceMin'  C     keep the position of this problem but do the stop later
234               ii = i
235               jj = j
236               icount = icount + 1
237            ENDIF            ENDIF
238    
239  C--   Fractional snow cover:  C--   Fractional snow cover:
# Line 267  C--   Compute ice temperatures Line 272  C--   Compute ice temperatures
272            tIc1(i,j) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1            tIc1(i,j) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
273            tIc2(i,j) = (Lfresh-qIc2(i,j)) / cpIce            tIc2(i,j) = (Lfresh-qIc2(i,j)) / cpIce
274    
275    #ifdef ALLOW_DBUG_THSICE
276            IF (tIc1(i,j).GT.0. _d 0 ) THEN            IF (tIc1(i,j).GT.0. _d 0 ) THEN
277             WRITE (standardMessageUnit,'(A,I12,1PE14.6)')             WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
278       &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)       &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)
# Line 279  C--   Compute ice temperatures Line 285  C--   Compute ice temperatures
285             WRITE (standardMessageUnit,'(A,4I5,2F11.4)')             WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
286       &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)       &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
287            ENDIF            ENDIF
 #ifdef ALLOW_DBUG_THSICE  
288            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
289       &      'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf(i,j),tIc1(i,j),tIc2(i,j)       &      'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf(i,j),tIc1(i,j),tIc2(i,j)
290  #endif  #endif
# Line 303  C--   Compute coefficients used in quadr Line 308  C--   Compute coefficients used in quadr
308          ENDIF          ENDIF
309         ENDDO         ENDDO
310        ENDDO        ENDDO
311          IF ( icount .gt. 0 ) THEN
312           WRITE (standardMessageUnit,'(A,I5,A)')
313         &      'THSICE_SOLVE4TEMP: there are ',icount,
314         &      ' case(s) where hIce<hIceMin;'
315           WRITE (standardMessageUnit,'(A,I3,A1,I3,A)')
316         &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
317         &      ') with hIce = ', hIce(ii,jj)
318           STOP 'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
319          ENDIF
320    
321  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
322  C--   Get surface fluxes over melting surface  C--   Get surface fluxes over melting surface
# Line 442  C     If no convergence, then repeat. Line 456  C     If no convergence, then repeat.
456  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
457            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
458       &    'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)       &    'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
 #endif  
459            IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN            IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN
460              WRITE (6,'(A,4I4,I12,F15.9)')             WRITE (6,'(A,4I4,I12,F15.9)')
461       &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)       &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
462              WRITE (6,*) 'BB: not converge: Tsf, dTsf=',             WRITE (6,*) 'BB: not converge: Tsf, dTsf=',
463       &                  Tsf(i,j), dTsrf(i,j)       &          Tsf(i,j), dTsrf(i,j)
464              WRITE (6,*) 'BB: not converge: flxNet,dFlxT=',             WRITE (6,*) 'BB: not converge: flxNet,dFlxT=',
465       &                  flxNet, dFlxdT(i,j)       &          flxNet, dFlxdT(i,j)
466              IF (Tsf(i,j).LT.-70. _d 0) STOP             IF (Tsf(i,j).LT.-70. _d 0) STOP
467            ENDIF            ENDIF
468    #endif
469    
470           ENDIF           ENDIF
471          ENDDO          ENDDO
472         ENDDO         ENDDO
   
473  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
474         ENDIF         ENDIF
475  #endif  #endif
# Line 511  C--   Compute new enthalpy for each laye Line 524  C--   Compute new enthalpy for each laye
524       &                + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))       &                + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))
525            qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh            qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh
526    
527    #ifdef ALLOW_DBUG_THSICE
528  C--   Make sure internal ice temperatures do not exceed Tmlt.  C--   Make sure internal ice temperatures do not exceed Tmlt.
529  C     (This should not happen for reasonable values of i0swFrac)  C     (This should not happen for reasonable values of i0swFrac)
530            IF (tIc1(i,j) .GE. Tmlt1) THEN            IF (tIc1(i,j) .GE. Tmlt1) THEN
# Line 522  C     (This should not happen for reason Line 536  C     (This should not happen for reason
536       &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)       &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)
537            ENDIF            ENDIF
538    
 #ifdef ALLOW_DBUG_THSICE  
539            IF ( dBug(i,j,bi,bj) ) THEN            IF ( dBug(i,j,bi,bj) ) THEN
540             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
541       &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)       &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)
Legend:
Removed from v.1.21  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22