/[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.24 by heimbach, Sat Oct 16 19:22:34 2010 UTC revision 1.27 by jmc, Fri Oct 22 13:43:35 2010 UTC
# Line 142  C     dt           :: timestep Line 142  C     dt           :: timestep
142        LOGICAL useBlkFlx        LOGICAL useBlkFlx
143        LOGICAL iterate4Tsf        LOGICAL iterate4Tsf
144        LOGICAL iceFlag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        LOGICAL iceFlag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145        INTEGER i,j        INTEGER i, j, k, iterMax
146        INTEGER k, iterMax, ii, jj, icount        INTEGER ii, jj, icount, stdUnit
147          CHARACTER*(MAX_LEN_MBUF) msgBuf
148        _RL  frsnow        _RL  frsnow
149        _RL  fswpen        _RL  fswpen
150        _RL  fswdn        _RL  fswdn
# Line 179  C-    Define grid-point location where t Line 180  C-    Define grid-point location where t
180  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
181    
182  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
183  c     act1 = bi - myBxLo(myThid)        act1 = bi - myBxLo(myThid)
184  c     max1 = myBxHi(myThid) - myBxLo(myThid) + 1        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
185  c     act2 = bj - myByLo(myThid)        act2 = bj - myByLo(myThid)
186  c     max2 = myByHi(myThid) - myByLo(myThid) + 1        max2 = myByHi(myThid) - myByLo(myThid) + 1
187  c     act3 = myThid - 1        act3 = myThid - 1
188  c     max3 = nTx*nTy        max3 = nTx*nTy
189  c     act4 = ikey_dynamics - 1        act4 = ikey_dynamics - 1
190        iicekey = (act1 + 1) + act2*max1        iicekey = (act1 + 1) + act2*max1
191       &                     + act3*max1*max2       &                     + act3*max1*max2
192       &                     + act4*max1*max2*max3       &                     + act4*max1*max2*max3
# Line 193  c     act4 = ikey_dynamics - 1 Line 194  c     act4 = ikey_dynamics - 1
194    
195  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
196  CADJ STORE flxsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE flxsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
197        DO j = jMin, jMax        DO j = 1-OLy, sNy+OLy
198         DO i = iMin, iMax         DO i = 1-OLx, sNx+OLx
199          tic1(i,j) = 0. _d 0          tIc1(i,j) = 0. _d 0
200          tic2(i,j) = 0. _d 0          tIc2(i,j) = 0. _d 0
201         END DO         ENDDO
202        END DO        ENDDO
203  #endif  #endif
204    
205          stdUnit = standardMessageUnit
206        useBlkFlx = useEXF .OR. useBulkForce        useBlkFlx = useEXF .OR. useBulkForce
207        dt  = thSIce_dtTemp        dt  = thSIce_dtTemp
208        IF ( dhSnowLin.GT.0. _d 0 ) THEN        IF ( dhSnowLin.GT.0. _d 0 ) THEN
# Line 211  CADJ STORE flxsw(:,:) = comlev1_bibj,key Line 213  CADJ STORE flxsw(:,:) = comlev1_bibj,key
213    
214        iterate4Tsf = .FALSE.        iterate4Tsf = .FALSE.
215        icount = 0        icount = 0
216  C      
217        DO j = jMin, jMax        DO j = jMin, jMax
218         DO i = iMin, iMax         DO i = iMin, iMax
219  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 223  c    &         + sNx*sNy*max1*max2*act3 Line 225  c    &         + sNx*sNy*max1*max2*act3
225  c    &         + sNx*sNy*max1*max2*max3*act4  c    &         + sNx*sNy*max1*max2*max3*act4
226  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
227  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
228  cCADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  devdt(i,j)        = comlev1_thsice_1, key=ikey_1
229  cCADJ STORE  dFlxdT        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  dFlxdT(i,j)        = comlev1_thsice_1, key=ikey_1
 cCADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1  
230  cCADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1  cCADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1
231  cCADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1
232  cCADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1
# Line 235  cCADJ STORE  tsrf(i,j)    = comlev1_thsi Line 236  cCADJ STORE  tsrf(i,j)    = comlev1_thsi
236            iterate4Tsf  = .TRUE.            iterate4Tsf  = .TRUE.
237            iceFlag(i,j) = .TRUE.            iceFlag(i,j) = .TRUE.
238  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
239            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,2I4,2I2)')
240       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
241  #endif  #endif
242            IF ( hIce(i,j).LT.hIceMin ) THEN            IF ( hIce(i,j).LT.hIceMin ) THEN
# Line 284  C--   Compute ice temperatures Line 285  C--   Compute ice temperatures
285    
286  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
287            IF (tIc1(i,j).GT.0. _d 0 ) THEN            IF (tIc1(i,j).GT.0. _d 0 ) THEN
288             WRITE (standardMessageUnit,'(A,I12,1PE14.6)')             WRITE(stdUnit,'(A,I12,1PE14.6)')
289       &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)       &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)
290             WRITE (standardMessageUnit,'(A,4I5,2F11.4)')             WRITE(stdUnit,'(A,4I5,2F11.4)')
291       &      ' 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)
292            ENDIF            ENDIF
293            IF ( tIc2(i,j).GT.0. _d 0) THEN            IF ( tIc2(i,j).GT.0. _d 0) THEN
294             WRITE (standardMessageUnit,'(A,I12,1PE14.6)')             WRITE(stdUnit,'(A,I12,1PE14.6)')
295       &       ' BBerr: Tice(2) > 0 ; it=', myIter, qIc2(i,j)       &       ' BBerr: Tice(2) > 0 ; it=', myIter, qIc2(i,j)
296             WRITE (standardMessageUnit,'(A,4I5,2F11.4)')             WRITE(stdUnit,'(A,4I5,2F11.4)')
297       &      ' 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)
298            ENDIF            ENDIF
299            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
300       &      '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)
301  #endif  #endif
302    
# Line 319  C--   Compute coefficients used in quadr Line 320  C--   Compute coefficients used in quadr
320         ENDDO         ENDDO
321        ENDDO        ENDDO
322        IF ( icount .gt. 0 ) THEN        IF ( icount .gt. 0 ) THEN
323         WRITE (standardMessageUnit,'(A,I5,A)')         WRITE(stdUnit,'(A,I5,A)')
324       &      'THSICE_SOLVE4TEMP: there are ',icount,       &      'THSICE_SOLVE4TEMP: there are ',icount,
325       &      ' case(s) where hIce<hIceMin;'       &      ' case(s) where hIce<hIceMin;'
326         WRITE (standardMessageUnit,'(A,I3,A1,I3,A)')         WRITE(stdUnit,'(A,I3,A1,I3,A)')
327       &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,       &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
328       &      ') with hIce = ', hIce(ii,jj)       &      ') with hIce = ', hIce(ii,jj)
329         STOP 'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'         WRITE( msgBuf, '(A)')
330         &      'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
331           CALL PRINT_ERROR( msgBuf , myThid )
332           STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
333        ENDIF        ENDIF
334    
335  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
336    
337  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
338  CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE devdt(:,:)    = comlev1_bibj,key=iicekey,byte=isbyte
339  CADJ STORE tsf      = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE tsf(:,:)      = comlev1_bibj,key=iicekey,byte=isbyte
340  #endif  #endif
341    
342  C--   Get surface fluxes over melting surface  C--   Get surface fluxes over melting surface
# Line 380  C     Tsfc converges. Line 384  C     Tsfc converges.
384        ENDIF        ENDIF
385    
386  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
387  CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE devdt(:,:)    = comlev1_bibj,key=iicekey,byte=isbyte
388  CADJ STORE dflxdt   = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE dflxdt(:,:)   = comlev1_bibj,key=iicekey,byte=isbyte
389  CADJ STORE flx0exsw = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE flx0exsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
390  CADJ STORE flxtexsw = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE flxtexsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
391  #endif  #endif
392    
393  C ----- begin iteration  -----  C ----- begin iteration  -----
# Line 423  C-    could add this "ifdef" to hide THS Line 427  C-    could add this "ifdef" to hide THS
427         ENDIF         ENDIF
428    
429  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
430  CADJ STORE devdt(i,j) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE devdt(:,:)    = comlev1_bibj,key=iicekey,byte=isbyte
431  CADJ STORE dflxdt(i,j) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE dflxdt(:,:)   = comlev1_bibj,key=iicekey,byte=isbyte
432  CADJ STORE flxtexsw(i,j) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE flxtexsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
433  CADJ STORE iceflag(i,j) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE iceflag(:,:)  = comlev1_bibj,key=iicekey,byte=isbyte
434  CADJ STORE tsf(i,j) = comlev1_bibj,key=iicekey,byte=isbyte  CADJ STORE tsf(:,:)      = comlev1_bibj,key=iicekey,byte=isbyte
435  #endif  #endif
436    
437  C--   Compute new top layer and surface temperatures.  C--   Compute new top layer and surface temperatures.
# Line 438  C     with different coefficients. Line 442  C     with different coefficients.
442           IF ( iceFlag(i,j) ) THEN           IF ( iceFlag(i,j) ) THEN
443            flxNet = sHeat(i,j) + flxTexSW(i,j)            flxNet = sHeat(i,j) + flxTexSW(i,j)
444  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
445            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
446       &     'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',       &     'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',
447       &      flxNet, dFlxdT(i,j), k12(i,j), k12(i,j)-dFlxdT(i,j)       &      flxNet, dFlxdT(i,j), k12(i,j), k12(i,j)-dFlxdT(i,j)
448  #endif  #endif
# Line 453  C     with different coefficients. Line 457  C     with different coefficients.
457            Tsf(i,j) = Tsf(i,j) + dTsrf(i,j)            Tsf(i,j) = Tsf(i,j) + dTsrf(i,j)
458            IF ( Tsf(i,j) .GT. 0. _d 0 ) THEN            IF ( Tsf(i,j) .GT. 0. _d 0 ) THEN
459  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
460             IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)             IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
461       &     '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)
462  #endif  #endif
463             a1 = a10(i,j) + k12(i,j)             a1 = a10(i,j) + k12(i,j)
# Line 482  C      if Tsfc has not converged.) Line 486  C      if Tsfc has not converged.)
486  C     If no convergence, then repeat.  C     If no convergence, then repeat.
487    
488  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
489            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
490       &    '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)
491            IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN            IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN
492             WRITE (6,'(A,4I4,I12,F15.9)')             WRITE(stdUnit,'(A,4I4,I12,F15.9)')
493       &       ' 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)
494             WRITE (6,*) 'BB: not converge: Tsf, dTsf=',             WRITE(stdUnit,*)
495       &          Tsf(i,j), dTsrf(i,j)       &        'BB: not converge: Tsf, dTsf=', Tsf(i,j), dTsrf(i,j)
496             WRITE (6,*) 'BB: not converge: flxNet,dFlxT=',             WRITE(stdUnit,*)
497       &          flxNet, dFlxdT(i,j)       &        'BB: not converge: flxNet,dFlxT=', flxNet, dFlxdT(i,j)
498             IF (Tsf(i,j).LT.-70. _d 0) STOP             IF ( Tsf(i,j).LT.-70. _d 0 ) THEN
499                 WRITE( msgBuf, '(A,2I4,2I3,I10,F12.3)')
500         &        'THSICE_SOLVE4TEMP: Too low Tsf in', i, j, bi, bj,
501         &                            myIter, Tsf(i,j)
502                 CALL PRINT_ERROR( msgBuf , myThid )
503                 STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
504               ENDIF
505            ENDIF            ENDIF
506  #endif  #endif
507    
# Line 516  C--   Compute new bottom layer temperatu Line 526  C--   Compute new bottom layer temperatu
526       &                 + rhoi*cpIce*hIce(i,j)*tIc2(i,j))       &                 + rhoi*cpIce*hIce(i,j)*tIc2(i,j))
527       &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))       &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))
528  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
529            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
530       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf(i,j),tIc1(i,j),tIc2(i,j)       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf(i,j),tIc1(i,j),tIc2(i,j)
531            netSW   = flxAtm(i,j)            netSW   = flxAtm(i,j)
532  #endif  #endif
# Line 542  C-    excess of energy @ surface (used f Line 552  C-    excess of energy @ surface (used f
552            sHeat(i,j) = flxNet - fct            sHeat(i,j) = flxNet - fct
553    
554  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
555            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
556       &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',       &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
557       &                    flxNet,fct,flxNet-fct,flxCnB(i,j)       &                    flxNet,fct,flxNet-fct,flxCnB(i,j)
558  #endif  #endif
# Line 556  C--   Compute new enthalpy for each laye Line 566  C--   Compute new enthalpy for each laye
566  C--   Make sure internal ice temperatures do not exceed Tmlt.  C--   Make sure internal ice temperatures do not exceed Tmlt.
567  C     (This should not happen for reasonable values of i0swFrac)  C     (This should not happen for reasonable values of i0swFrac)
568            IF (tIc1(i,j) .GE. Tmlt1) THEN            IF (tIc1(i,j) .GE. Tmlt1) THEN
569             WRITE (6,'(A,2I4,2I3,1P2E14.6)')             WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
570       &     ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,tIc1(i,j),Tmlt1       &     ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,tIc1(i,j),Tmlt1
571            ENDIF            ENDIF
572            IF (tIc2(i,j) .GE. 0. _d 0) THEN            IF (tIc2(i,j) .GE. 0. _d 0) THEN
573             WRITE (6,'(A,2I4,2I3,1P2E14.6)')             WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
574       &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)       &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)
575            ENDIF            ENDIF
576    
577            IF ( dBug(i,j,bi,bj) ) THEN            IF ( dBug(i,j,bi,bj) ) THEN
578             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',             WRITE(stdUnit,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
579       &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)       &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)
580             WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',             WRITE(stdUnit,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
581       &           sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)       &           sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)
582             WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',             WRITE(stdUnit,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
583       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
584            ENDIF            ENDIF
585  #endif  #endif
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22