/[MITgcm]/MITgcm/pkg/seaice/seaice_solve4temp.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_solve4temp.F

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

revision 1.1 by jmc, Thu Sep 23 22:46:24 2010 UTC revision 1.2 by jmc, Fri Sep 24 13:52:07 2010 UTC
# Line 106  c     Constants to calculate Saturation Line 106  c     Constants to calculate Saturation
106        _RL A3        (1:sNx,1:sNy)        _RL A3        (1:sNx,1:sNy)
107        _RL B         (1:sNx,1:sNy)        _RL B         (1:sNx,1:sNy)
108        _RL A1        (1:sNx,1:sNy)        _RL A1        (1:sNx,1:sNy)
109  #else  #else  /* USE_ORIGINAL_SBI */
110        _RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t        _RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t
111  C     specific humidity at ice surface variables  C     specific humidity at ice surface variables
112        _RL  mm_pi,mm_log10pi,dqhice_dTice        _RL  mm_pi,mm_log10pi,dqhice_dTice
113  #endif  #endif /* USE_ORIGINAL_SBI */
114    
115  C     effective conductivity of combined ice and snow  C     effective conductivity of combined ice and snow
116        _RL  effConduct        _RL  effConduct
# Line 132  C MAYKUTS CONSTANTS FOR SAT. VAP. PRESSU Line 132  C MAYKUTS CONSTANTS FOR SAT. VAP. PRESSU
132    
133        QS1=0.622 _d +00/1013.0 _d +00        QS1=0.622 _d +00/1013.0 _d +00
134    
135  #else  #else /* USE_ORIGINAL_SBI */
136        aa1 = 2663.5 _d 0        aa1 = 2663.5 _d 0
137        aa2 = 12.537 _d 0        aa2 = 12.537 _d 0
138        bb1 = 0.622 _d 0        bb1 = 0.622 _d 0
# Line 141  C MAYKUTS CONSTANTS FOR SAT. VAP. PRESSU Line 141  C MAYKUTS CONSTANTS FOR SAT. VAP. PRESSU
141        cc0 = TEN ** aa2        cc0 = TEN ** aa2
142        cc1 = cc0*aa1*bb1*Ppascals*log(10. _d 0)        cc1 = cc0*aa1*bb1*Ppascals*log(10. _d 0)
143        cc2 = cc0*bb2        cc2 = cc0*bb2
144  #endif  #endif /* USE_ORIGINAL_SBI */
145    
146  C     FREEZING TEMPERATURE OF SEAWATER  C     FREEZING TEMPERATURE OF SEAWATER
147  #ifndef SEAICE_VARIABLE_FREEZING_POINT  #ifdef SEAICE_VARIABLE_FREEZING_POINT
148    c     Use a variable seawater freezing point
149             TB = -0.0575 _d 0*salt(I,J,kSrf,bi,bj) + 0.0901 _d 0
150         &        + 273.15 _d 0
151    #else /* SEAICE_VARIABLE_FREEZING_POINT */
152  c     Use a constant seaswater freezing point  c     Use a constant seaswater freezing point
153  #ifdef USE_ORIGINAL_SBI  #ifdef USE_ORIGINAL_SBI
154        TB=271.2 _d 0        TB=271.2 _d 0
155  #else  #else /* USE_ORIGINAL_SBI */
156        TB=273.15 _d 0 + SEAICE_freeze        TB=273.15 _d 0 + SEAICE_freeze
157  #endif  #endif /* USE_ORIGINAL_SBI */
158  #else  #endif /* SEAICE_VARIABLE_FREEZING_POINT */
 c     Use a variable seawater freezing point  
          TB = -0.0575 _d 0*salt(I,J,kSrf,bi,bj) + 0.0901 _d 0  
      &        + 273.15 _d 0  
 #endif  
159    
160  C     SENSIBLE HEAT CONSTANT  C     SENSIBLE HEAT CONSTANT
161        D1=SEAICE_sensHeat        D1=SEAICE_sensHeat
# Line 171  C     MELTING TEMPERATURE OF ICE Line 171  C     MELTING TEMPERATURE OF ICE
171        TMELT=273.16 _d +00        TMELT=273.16 _d +00
172        TMELTP=273.159 _d +00        TMELTP=273.159 _d +00
173        SurfMeltTemp = TMELTP        SurfMeltTemp = TMELTP
174  #else  #else /* USE_ORIGINAL_SBI */
175        TMELT = 273.15 _d 0        TMELT = 273.15 _d 0
176        SurfMeltTemp = TMELT        SurfMeltTemp = TMELT
177  #endif  #endif /* USE_ORIGINAL_SBI */
178    
179  C     ICE CONDUCTIVITY  C     ICE CONDUCTIVITY
180        XKI=SEAICE_iceConduct        XKI=SEAICE_iceConduct
# Line 211  c           Reset the snow/ice surface t Line 211  c           Reset the snow/ice surface t
211              A3(I,J) = ZERO              A3(I,J) = ZERO
212  c            B(I,J) = ZERO  c            B(I,J) = ZERO
213              lwdownLoc(I,J) = MAX(MIN_LWDOWN,LWDOWN(I,J,bi,bj))              lwdownLoc(I,J) = MAX(MIN_LWDOWN,LWDOWN(I,J,bi,bj))
214  #else  #else /* USE_ORIGINAL_SBI */
215              F_swi    (I,J) = 0. _d 0              F_swi    (I,J) = 0. _d 0
216              F_lwd    (I,J) = 0. _d 0              F_lwd    (I,J) = 0. _d 0
217              F_lwu    (I,J) = 0. _d 0              F_lwu    (I,J) = 0. _d 0
# Line 221  c            B(I,J) = ZERO Line 221  c            B(I,J) = ZERO
221              tsurfLoc(I,J) =  TSURF(I,J,bi,bj)              tsurfLoc(I,J) =  TSURF(I,J,bi,bj)
222              atempLoc (I,J) = MAX(TMELT + MIN_ATEMP,ATEMP(I,J,bi,bj))              atempLoc (I,J) = MAX(TMELT + MIN_ATEMP,ATEMP(I,J,bi,bj))
223              lwdownLoc(I,J) = LWDOWN(I,J,bi,bj)              lwdownLoc(I,J) = LWDOWN(I,J,bi,bj)
224  #endif  #endif /* USE_ORIGINAL_SBI */
225    
226           ENDDO           ENDDO
227        ENDDO        ENDDO
# Line 264  c              (What is the source of th Line 264  c              (What is the source of th
264       &                            ALB_SNOW(I,J))       &                            ALB_SNOW(I,J))
265                 ENDIF                 ENDIF
266    
267  #else  #else /* USE_ORIGINAL_SBI */
268                 IF (HSNOW_ACTUAL(I,J) .GT. ZERO) THEN                 IF (HSNOW_ACTUAL(I,J) .GT. ZERO) THEN
269                    ALB(I,J) = ALB_SNOW(I,J)                    ALB(I,J) = ALB_SNOW(I,J)
270                 ELSE                 ELSE
271                    ALB(I,J) = ALB_ICE(I,J)                    ALB(I,J) = ALB_ICE(I,J)
272                 ENDIF                 ENDIF
273  #endif  #endif /* USE_ORIGINAL_SBI */
274    
275    
276  #ifdef USE_ORIGINAL_SBI  #ifdef USE_ORIGINAL_SBI
# Line 289  C        SW PENETRATION UNDER ICE Line 289  C        SW PENETRATION UNDER ICE
289       &        +lwdownLoc(I,J)*0.97 _d 0       &        +lwdownLoc(I,J)*0.97 _d 0
290       &        +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)       &        +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
291          ENDIF          ENDIF
292  #endif  #endif /* ALLOW_DOWNWARD_RADIATION */
293    
294  #else  #else /* USE_ORIGINAL_SBI */
295    
296  c     The longwave radiative flux convergence  c     The longwave radiative flux convergence
297                 F_lwd(I,J) = - 0.97 _d 0 * lwdownLoc(I,J)                 F_lwd(I,J) = - 0.97 _d 0 * lwdownLoc(I,J)
# Line 325  c     Set a mininum sea ice thickness of Line 325  c     Set a mininum sea ice thickness of
325  c     the magnitude of conductive heat fluxes.  c     the magnitude of conductive heat fluxes.
326                 HICE_ACTUAL(I,J) = max(HICE_ACTUAL(I,J),5. _d -2)                 HICE_ACTUAL(I,J) = max(HICE_ACTUAL(I,J),5. _d -2)
327    
328  #endif  #endif /* USE_ORIGINAL_SBI */
329    
330  c     The effective conductivity of the two-layer  c     The effective conductivity of the two-layer
331  c     snow/ice system.  c     snow/ice system.
# Line 333  c     snow/ice system. Line 333  c     snow/ice system.
333                 effConduct=                 effConduct=
334       &             XKS/(HSNOW_ACTUAL(I,J)/HICE_ACTUAL(I,J) +       &             XKS/(HSNOW_ACTUAL(I,J)/HICE_ACTUAL(I,J) +
335       &                  XKS/XKI)/HICE_ACTUAL(I,J)       &                  XKS/XKI)/HICE_ACTUAL(I,J)
336  #else  #else /* USE_ORIGINAL_SBI */
337                 effConduct = XKI * XKS /                 effConduct = XKI * XKS /
338       &            (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))       &            (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))
339  #endif  #endif /* USE_ORIGINAL_SBI */
340    
341    
342    
# Line 366  c     snow/ice system. Line 366  c     snow/ice system.
366                    print '(A,i6)','ibi energy balance iterat ', myIter                    print '(A,i6)','ibi energy balance iterat ', myIter
367    
368                 ENDIF                 ENDIF
369  #endif  #endif /* SEAICE_DEBUG */
370    
371  ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
372                 DO ITER=1,IMAX_TICE                 DO ITER=1,IMAX_TICE
# Line 381  c                 Calculate the specific Line 381  c                 Calculate the specific
381  c                 Use the Maykut polynomial  c                 Use the Maykut polynomial
382                    qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)                    qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
383    
384  #else  #else /* USE_ORIGINAL_SBI */
385  c                 Use an approximation which is more accurate at low temperatures  c                 Use an approximation which is more accurate at low temperatures
386    
387  c                 log 10 of the sat vap pressure  c                 log 10 of the sat vap pressure
# Line 393  c                 boundary layer (BL) ab Line 393  c                 boundary layer (BL) ab
393    
394                    qhice(I,J) = bb1*mm_pi / (Ppascals - (ONE - bb1) *                    qhice(I,J) = bb1*mm_pi / (Ppascals - (ONE - bb1) *
395       &               mm_pi)       &               mm_pi)
396  #endif  #endif /* USE_ORIGINAL_SBI */
397    
398  c                 Caclulate the flux terms based on the updated tsurfLoc  c                 Caclulate the flux terms based on the updated tsurfLoc
399  #ifdef USE_ORIGINAL_SBI  #ifdef USE_ORIGINAL_SBI
400                    A2(I,J)=-D1*UG(I,J)*t1-D1I*UG(I,J)*qhice(I,J)-D3*t4                    A2(I,J)=-D1*UG(I,J)*t1-D1I*UG(I,J)*qhice(I,J)-D3*t4
401                    A3(I,J) = 4.0 _d 0 * D3 * t3 + effConduct + D1*UG(I,J)                    A3(I,J) = 4.0 _d 0 * D3 * t3 + effConduct + D1*UG(I,J)
402                    F_c(I,J)=-effConduct*(TB-tsurfLoc(I,J))                    F_c(I,J)=-effConduct*(TB-tsurfLoc(I,J))
403  #else  #else /* USE_ORIGINAL_SBI */
404  c                 A constant for SVP derivative w.r.t TICE  c                 A constant for SVP derivative w.r.t TICE
405                    cc3t = TEN **(aa1 / t1)                    cc3t = TEN **(aa1 / t1)
406    
# Line 422  c                 d(F_ia)/d(TICE) Line 422  c                 d(F_ia)/d(TICE)
422                    F_ia(I,J)  = F_lwd(I,J) + F_swi(I,J) + F_lwu(I,J) +                    F_ia(I,J)  = F_lwd(I,J) + F_swi(I,J) + F_lwu(I,J) +
423       &               F_c(I,J) + F_sens(I,J) + F_lh(I,J)       &               F_c(I,J) + F_sens(I,J) + F_lh(I,J)
424    
425  #endif  #endif /* USE_ORIGINAL_SBI */
426    
427  #ifdef SEAICE_DEBUG  #ifdef SEAICE_DEBUG
428                    IF ( (I .EQ. SEAICE_debugPointX)   .and.                    IF ( (I .EQ. SEAICE_debugPointX)   .and.
# Line 438  c                 d(F_ia)/d(TICE) Line 438  c                 d(F_ia)/d(TICE)
438                       print '(A,i6,4(1x,D24.15))',                       print '(A,i6,4(1x,D24.15))',
439       &                  'ice-iter A3 (-A1+A2)  ', ITER, A3(I,J),       &                  'ice-iter A3 (-A1+A2)  ', ITER, A3(I,J),
440       &                  -(A1(I,J) + A2(I,J))       &                  -(A1(I,J) + A2(I,J))
441  #else  #else /* USE_ORIGINAL_SBI */
442    
443                       print '(A,i6,4(1x,D24.15))',                       print '(A,i6,4(1x,D24.15))',
444       &                  'ice-iter dFiDTs1 F_ia ', ITER, dFiDTs1,       &                  'ice-iter dFiDTs1 F_ia ', ITER, dFiDTs1,
445       &                  F_ia(I,J)       &                  F_ia(I,J)
446  #endif  #endif /* USE_ORIGINAL_SBI */
447    
448                    ENDIF                    ENDIF
449  #endif  #endif /* SEAICE_DEBUG */
450    
451  c                 Update tsurfLoc  c                 Update tsurfLoc
452  #ifdef USE_ORIGINAL_SBI  #ifdef USE_ORIGINAL_SBI
# Line 456  c                 Update tsurfLoc Line 456  c                 Update tsurfLoc
456                    tsurfLoc(I,J) =MAX(273.16 _d 0+MIN_TICE,tsurfLoc(I,J))                    tsurfLoc(I,J) =MAX(273.16 _d 0+MIN_TICE,tsurfLoc(I,J))
457                    tsurfLoc(I,J) =MIN(tsurfLoc(I,J),TMELT)                    tsurfLoc(I,J) =MIN(tsurfLoc(I,J),TMELT)
458    
459  #else  #else /* USE_ORIGINAL_SBI */
460                    tsurfLoc(I,J) = tsurfLoc(I,J) - F_ia(I,J) / dFiDTs1                    tsurfLoc(I,J) = tsurfLoc(I,J) - F_ia(I,J) / dFiDTs1
461    
462  c                 If the search leads to tsurfLoc < 50 Kelvin,  c                 If the search leads to tsurfLoc < 50 Kelvin,
# Line 468  c                 realistic values. Line 468  c                 realistic values.
468                    IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) THEN                    IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) THEN
469                       tsurfLoc(I,J) = TMELT                       tsurfLoc(I,J) = TMELT
470                    ENDIF                    ENDIF
471  #endif  #endif /* USE_ORIGINAL_SBI */
472    
473  #ifdef SEAICE_DEBUG  #ifdef SEAICE_DEBUG
474                    IF ( (I .EQ. SEAICE_debugPointX)   .and.                    IF ( (I .EQ. SEAICE_debugPointX)   .and.
# Line 479  c                 realistic values. Line 479  c                 realistic values.
479       &                  tsurfLoc(I,J),       &                  tsurfLoc(I,J),
480       &                  log10(abs(tsurfLoc(I,J) - t1))       &                  log10(abs(tsurfLoc(I,J) - t1))
481                    ENDIF                    ENDIF
482  #endif  #endif /* SEAICE_DEBUG */
483    
484                 ENDDO            !/* Iterations */                 ENDDO            !/* Iterations */
485  ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
# Line 498  C              SW PENETRATION UNDER ICE Line 498  C              SW PENETRATION UNDER ICE
498  #ifdef ALLOW_DOWNWARD_RADIATION  #ifdef ALLOW_DOWNWARD_RADIATION
499                     IcePenetSWFlux(I,J)=-(ONE-ALB(I,J))*SWDOWN(I,J,bi,bj)                     IcePenetSWFlux(I,J)=-(ONE-ALB(I,J))*SWDOWN(I,J,bi,bj)
500       &              *XIO*EXP(-1.5 _d 0*HICE_ACTUAL(I,J))       &              *XIO*EXP(-1.5 _d 0*HICE_ACTUAL(I,J))
501  #endif  #endif /* ALLOW_DOWNWARD_RADIATION */
502                 ENDIF                 ENDIF
503    
504  #else  #else /* USE_ORIGINAL_SBI */
505                 tsurfLoc(I,J) = MIN(tsurfLoc(I,J),TMELT)                 tsurfLoc(I,J) = MIN(tsurfLoc(I,J),TMELT)
506                 TSURF(I,J,bi,bj) = tsurfLoc(I,J)                 TSURF(I,J,bi,bj) = tsurfLoc(I,J)
507    
# Line 529  c              The flux between the ice/ Line 529  c              The flux between the ice/
529  c              (excludes upward conductive fluxes)  c              (excludes upward conductive fluxes)
530                 F_ia(I,J)    = F_lwd(I,J) + F_swi(I,J) + F_lwu(I,J) +                 F_ia(I,J)    = F_lwd(I,J) + F_swi(I,J) + F_lwu(I,J) +
531       &            F_sens(I,J) + F_lh(I,J)       &            F_sens(I,J) + F_lh(I,J)
532  #endif  #endif /* USE_ORIGINAL_SBI */
533    
534  c              Caclulate the net ice-ocean and ice-atmosphere fluxes  c              Caclulate the net ice-ocean and ice-atmosphere fluxes
535                 IF (F_c(I,J) .LT. ZERO) THEN                 IF (F_c(I,J) .LT. ZERO) THEN
# Line 537  c              Caclulate the net ice-oce Line 537  c              Caclulate the net ice-oce
537                    F_ia_net(I,J) = ZERO                    F_ia_net(I,J) = ZERO
538                 ELSE                 ELSE
539                    F_io_net(I,J) = ZERO                    F_io_net(I,J) = ZERO
540                    F_ia_net(I,J) = F_lwd(I,J) + F_swi(I,J) + F_lwu(I,J) +                    F_ia_net(I,J) = F_ia(I,J)
      &               F_sens(I,J) + F_lh(I,J)  
541                 ENDIF  !/* conductive fluxes up or down */                 ENDIF  !/* conductive fluxes up or down */
542    
543    
# Line 578  c     &               F_c(I,J), F_lh(I,J Line 577  c     &               F_c(I,J), F_lh(I,J
577       &                -(A1(I,J)+A2(I,J)),       &                -(A1(I,J)+A2(I,J)),
578       &                -(A1(I,J)+A2(I,J)-F_c(I,J)),       &                -(A1(I,J)+A2(I,J)-F_c(I,J)),
579       &                F_c(I,J)       &                F_c(I,J)
580  #else  #else /* USE_ORIGINAL_SBI */
581       &                F_ia(I,J),       &                F_ia(I,J),
582       &                F_ia_net(I,J),       &                F_ia_net(I,J),
583       &                F_c(I,J)       &                F_c(I,J)
584  #endif  #endif /* USE_ORIGINAL_SBI */
585    
586                    print '(A)','----------------------------------------'                    print '(A)','----------------------------------------'
587    
588                 ENDIF                 ENDIF
589  #endif  #endif /* SEAICE_DEBUG */
590    
591              ENDIF               !/* HICE_ACTUAL > 0 */              ENDIF               !/* HICE_ACTUAL > 0 */
592    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22