--- MITgcm/pkg/aim_v23/phy_suflux_sice.F	2004/05/21 17:43:04	1.3
+++ MITgcm/pkg/aim_v23/phy_suflux_sice.F	2004/06/24 23:43:11	1.4
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/aim_v23/phy_suflux_sice.F,v 1.3 2004/05/21 17:43:04 jmc Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/aim_v23/phy_suflux_sice.F,v 1.4 2004/06/24 23:43:11 jmc Exp $
 C $Name:  $
 
 #include "AIM_OPTIONS.h"
@@ -9,9 +9,9 @@
       SUBROUTINE SUFLUX_SICE(
      I                   PSA, FMASK, EMISloc,
      I                   Tsurf, dTskin, SSR, SLRD,
-     I                   T0, Q0, CDENVV,
+     I                   T1, T0, Q0, DENVV,
      O                   SHF, EVAP, SLRU,
-     O                   Evp0, dEvp, Slr0, dSlr, sFlx,
+     O                   Shf0, dShf, Evp0, dEvp, Slr0, dSlr, sFlx,
      O                   TSFC, TSKIN,
      I                   bi,bj,myThid)
 
@@ -52,13 +52,16 @@
 C    dTskin :: temp. correction for daily-cycle heating [K]
 C    SSR    :: sfc sw radiation (net flux)     (2-dim)
 C    SLRD   :: sfc lw radiation (downward flux)(2-dim)
+C    T1     :: near-surface air temperature (from Pot.temp)
 C    T0     :: near-surface air temperature    (2-dim)
 C    Q0     :: near-surface sp. humidity [g/kg](2-dim)
-C    CDENVV :: sensible heat flux coefficient  (2-dim)
+C    DENVV  :: surface flux (sens,lat.) coeff. (=Rho*|V|) [kg/m2/s]
 C--   Output:
 C    SHF    :: sensible heat flux              (2-dim)
 C    EVAP   :: evaporation [g/(m^2 s)]         (2-dim)
 C    SLRU   :: sfc lw radiation (upward flux)  (2-dim)
+C    Shf0   :: sensible heat flux over freezing surf.
+C    dShf   :: sensible heat flux derivative relative to surf. temp
 C    Evp0   :: evaporation computed over freezing surface (Ts=0.oC)
 C    dEvp   :: evaporation derivative relative to surf. temp
 C    Slr0   :: upward long wave radiation over freezing surf.
@@ -74,10 +77,11 @@
       _RL  PSA(NGP), FMASK(NGP), EMISloc
       _RL  Tsurf(NGP), dTskin(NGP)
       _RL  SSR(NGP), SLRD(NGP)
-      _RL  T0(NGP), Q0(NGP), CDENVV(NGP)
+      _RL  T1(NGP), T0(NGP), Q0(NGP), DENVV(NGP)
 
       _RL  SHF(NGP), EVAP(NGP), SLRU(NGP)
-      _RL  Evp0(NGP), dEvp(NGP), Slr0(NGP), dSlr(NGP), sFlx(NGP,0:2)
+      _RL  Shf0(NGP), dShf(NGP), Evp0(NGP), dEvp(NGP)
+      _RL  Slr0(NGP), dSlr(NGP), sFlx(NGP,0:2)
       _RL  TSFC(NGP), TSKIN(NGP)
 
       INTEGER bi,bj,myThid
@@ -86,6 +90,9 @@
 #ifdef ALLOW_AIM
 
 C-- Local variables:
+C    CDENVV :: surf. heat flux (sens.,lat.) coeff including stability effect
+      _RL CDENVV(NGP), RDTH, FSSICE
+      _RL Fstb0, dTstb, dFstb
       _RL QSAT0(NGP,2)
       _RL QDUMMY(1), RDUMMY(1), TS2
       INTEGER J
@@ -105,29 +112,69 @@
 
 C--   2. Computation of fluxes over land and sea
 
-C     2.1 Wind stress
+C     2.1 Stability correction
 
-C     2.2 Sensible heat flux (from clim. TS over land)
+      RDTH = FSTAB/DTHETA
 
       DO J=1,NGP
-        SHF(J)   =  CDENVV(J)*CP*(TSKIN(J)-T0(J))
-        sFlx(J,0)= -CDENVV(J)*CP*(TSFC(J) -T0(J))
-        sFlx(J,1)= -SHF(J)
-        sFlx(J,2)= -CDENVV(J)*CP
+        FSSICE=1.+MIN(DTHETA,MAX(-DTHETA,TSKIN(J)-T1(J)))*RDTH
+        CDENVV(J)=CHS*DENVV(J)*FSSICE
       ENDDO
 
-C     2.3 Evaporation
+      IF ( dTstab.GT.0. _d 0 ) THEN
+C-    account for stability function derivative relative to Tsurf:
+C note: to avoid discontinuity in the derivative (because of min,max), compute
+C   the derivative using the discrete form: F(Ts+dTstab)-F(Ts-dTstab)/2.dTstab
+       DO J=1,NGP
+        Fstb0 = 1.+MIN(DTHETA,MAX(-DTHETA,TSFC(J) -T1(J)))*RDTH
+        Shf0(J) = CHL*DENVV(J)*Fstb0
+        dTstb = ( DTHETA+dTstab-ABS(TSKIN(J)-T1(J)) )/dTstab
+        dFstb = RDTH*MIN(1. _d 0, MAX(0. _d 0, dTstb*0.5 _d 0))
+        dShf(J) = CHL*DENVV(J)*dFstb
+       ENDDO
+      ENDIF
+
+C     2.2 Evaporation
 
       CALL SHTORH (2, NGP, TSKIN, PSA, 1. _d 0, QDUMMY, dEvp,
      &                QSAT0(1,1), myThid)
       CALL SHTORH (0, NGP, TSFC, PSA, 1. _d 0, QDUMMY, RDUMMY,
      &                QSAT0(1,2), myThid)
 
-      DO J=1,NGP
+      IF ( dTstab.GT.0. _d 0 ) THEN
+C-    account for stability function derivative relative to Tsurf:
+       DO J=1,NGP
+        EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
+        Evp0(J) =   Shf0(J)*(QSAT0(J,2)-Q0(J))
+        dEvp(J) = CDENVV(J)*dEvp(J)
+     &            + dShf(J)*(QSAT0(J,1)-Q0(J))
+       ENDDO
+      ELSE
+       DO J=1,NGP
         EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
         Evp0(J) = CDENVV(J)*(QSAT0(J,2)-Q0(J))
         dEvp(J) = CDENVV(J)*dEvp(J)
-      ENDDO
+       ENDDO
+      ENDIF
+
+C     2.3 Sensible heat flux
+
+      IF ( dTstab.GT.0. _d 0 ) THEN
+C-    account for stability function derivative relative to Tsurf:
+       DO J=1,NGP
+        SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
+        Shf0(J) =   Shf0(J)*CP*(TSFC(J) -T0(J))
+        dShf(J) = CDENVV(J)*CP
+     &            + dShf(J)*CP*(TSKIN(J)-T0(J))
+       ENDDO
+      ELSE
+       DO J=1,NGP
+        SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
+        Shf0(J) = CDENVV(J)*CP*(TSFC(J) -T0(J))
+        dShf(J) = CDENVV(J)*CP
+        dShf(J) = MAX( dShf(J), 0. _d 0 )
+       ENDDO
+      ENDIF
 
 C     2.4 Emission of lw radiation from the surface
 
@@ -141,12 +188,12 @@
 
 C--   Compute net surface heat flux and its derivative ./. surf. temp.
       DO J=1,NGP
-        sFlx(J,0)= sFlx(J,0)
-     &           - ALHC*Evp0(J) - EMISloc*Slr0(J) + SLRD(J)
-        sFlx(J,1)= sFlx(J,1)
-     &           - ALHC*EVAP(J) - EMISloc*SLRU(J) + SLRD(J)
-        sFlx(J,2)= sFlx(J,2)
-     &           - ALHC*dEvp(J) - EMISloc*dSlr(J)
+        sFlx(J,0)= ( SLRD(J) - EMISloc*Slr0(J) )
+     &           - ( Shf0(J) + ALHC*Evp0(J) )
+        sFlx(J,1)= ( SLRD(J) - EMISloc*SLRU(J) )
+     &           - ( SHF(J) + ALHC*EVAP(J) )
+        sFlx(J,2)=            -EMISloc*dSlr(J)
+     &           - ( dShf(J) + ALHC*dEvp(J) )
       ENDDO
       IF ( aim_energPrecip ) THEN
 C-     Evap of snow/ice: substract Latent Heat of freezing from heatFlux