pkg/aim_v23/phy_suflux_sice.F

C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/phy_suflux_sice.F,v 1.3 2004/05/21 17:43:04 jmc Exp $
C $Name:  $

#include "AIM_OPTIONS.h"

CBOP
C     !ROUTINE: SUFLUX_SICE
C     !INTERFACE:
      SUBROUTINE SUFLUX_SICE(
     I                   PSA, FMASK, EMISloc,
     I                   Tsurf, dTskin, SSR, SLRD,
     I                   T1, T0, Q0, DENVV,
     O                   SHF, EVAP, SLRU,
     O                   Shf0, dShf, Evp0, dEvp, Slr0, dSlr, sFlx,
     O                   TSFC, TSKIN,
     I                   bi,bj,myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R SUFLUX_SICE
C     | o compute surface flux over sea-ice
C     *==========================================================*
C     | o contains part of original S/R SUFLUX (Speedy code)
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     Resolution parameters

C-- size for MITgcm & Physics package :
#include "AIM_SIZE.h"
#include "EEPARAMS.h"

C-- Physics package
#include "AIM_PARAMS.h"

C     Physical constants + functions of sigma and latitude
#include "com_physcon.h"

C     Surface flux constants
#include "com_sflcon.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C--   Input:
C    PSA    :: norm. surface pressure [p/p0]   (2-dim)
C    FMASK  :: fractional land-sea mask        (2-dim)
C    EMISloc:: longwave surface emissivity
C    Tsurf  :: surface temperature        (2-dim)
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    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.
C    dSlr   :: upward long wave rad. derivative relative to surf. temp
C    sFlx   :: net heat flux (+=down) except SW, function of surf. temp Ts:
C              0: Flux(Ts=0.oC) ; 1: Flux(Ts^n) ; 2: d.Flux/d.Ts(Ts^n)
C    TSFC   :: surface temperature (clim.)     (2-dim)
C    TSKIN  :: skin surface temperature        (2-dim)
C--   Input:
C    bi,bj  :: tile index
C    myThid :: Thread number for this instance of the routine
C--
      _RL  PSA(NGP), FMASK(NGP), EMISloc
      _RL  Tsurf(NGP), dTskin(NGP)
      _RL  SSR(NGP), SLRD(NGP)
      _RL  T1(NGP), T0(NGP), Q0(NGP), DENVV(NGP)

      _RL  SHF(NGP), EVAP(NGP), SLRU(NGP)
      _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
CEOP

#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

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C     1.5 Define effective skin temperature to compensate for
C         non-linearity of heat/moisture fluxes during the daily cycle

      DO J=1,NGP
c       TSKIN(J) = Tsurf(J) + dTskin(J)
c       TSFC(J)=273.16 _d 0 + dTskin(J)
        TSKIN(J) = Tsurf(J)
        TSFC(J)=273.16 _d 0
      ENDDO


C--   2. Computation of fluxes over land and sea

C     2.1 Stability correction

      RDTH = FSTAB/DTHETA

      DO J=1,NGP
        FSSICE=1.+MIN(DTHETA,MAX(-DTHETA,TSKIN(J)-T1(J)))*RDTH
        CDENVV(J)=CHS*DENVV(J)*FSSICE
      ENDDO

      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)

      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
      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

      DO J=1,NGP
        TS2     = TSFC(J)*TSFC(J)
        Slr0(J) = SBC*TS2*TS2
        TS2     = TSKIN(J)*TSKIN(J)
        SLRU(J) = SBC*TS2*TS2
        dSlr(J)  = 4. _d 0 *SBC*TS2*TSKIN(J)
      ENDDO

C--   Compute net surface heat flux and its derivative ./. surf. temp.
      DO J=1,NGP
        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
       DO J=1,NGP
         sFlx(J,0) = sFlx(J,0) - ALHF*Evp0(J)
         sFlx(J,1) = sFlx(J,1) - ALHF*EVAP(J)
         sFlx(J,2) = sFlx(J,2) - ALHF*dEvp(J)
       ENDDO
      ENDIF

C--   3. Adjustment of skin temperature and fluxes over land
C--      based on energy balance (to be implemented)
C        <= done separately for each surface type (land,ocean,sea-ice)

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* ALLOW_AIM */

      RETURN
      END
1	jmc	1.4	C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/phy_suflux_sice.F,v 1.3 2004/05/21 17:43:04 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "AIM_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: SUFLUX_SICE
8			C !INTERFACE:
9			SUBROUTINE SUFLUX_SICE(
10			I PSA, FMASK, EMISloc,
11			I Tsurf, dTskin, SSR, SLRD,
12	jmc	1.4	I T1, T0, Q0, DENVV,
13	jmc	1.1	O SHF, EVAP, SLRU,
14	jmc	1.4	O Shf0, dShf, Evp0, dEvp, Slr0, dSlr, sFlx,
15	jmc	1.1	O TSFC, TSKIN,
16			I bi,bj,myThid)
17
18			C !DESCRIPTION: \bv
19			C ==========================================================
20			C \| S/R SUFLUX_SICE
21			C \| o compute surface flux over sea-ice
22			C ==========================================================
23			C \| o contains part of original S/R SUFLUX (Speedy code)
24			C ==========================================================
25			C \ev
26
27			C !USES:
28			IMPLICIT NONE
29
30			C Resolution parameters
31
32			C-- size for MITgcm & Physics package :
33			#include "AIM_SIZE.h"
34			#include "EEPARAMS.h"
35
36			C-- Physics package
37			#include "AIM_PARAMS.h"
38
39			C Physical constants + functions of sigma and latitude
40			#include "com_physcon.h"
41
42			C Surface flux constants
43			#include "com_sflcon.h"
44
45			C !INPUT/OUTPUT PARAMETERS:
46			C == Routine Arguments ==
47			C-- Input:
48			C PSA :: norm. surface pressure [p/p0] (2-dim)
49			C FMASK :: fractional land-sea mask (2-dim)
50			C EMISloc:: longwave surface emissivity
51			C Tsurf :: surface temperature (2-dim)
52			C dTskin :: temp. correction for daily-cycle heating [K]
53			C SSR :: sfc sw radiation (net flux) (2-dim)
54			C SLRD :: sfc lw radiation (downward flux)(2-dim)
55	jmc	1.4	C T1 :: near-surface air temperature (from Pot.temp)
56	jmc	1.1	C T0 :: near-surface air temperature (2-dim)
57			C Q0 :: near-surface sp. humidity [g/kg](2-dim)
58	jmc	1.4	C DENVV :: surface flux (sens,lat.) coeff. (=Rho*\|V\|) [kg/m2/s]
59	jmc	1.1	C-- Output:
60			C SHF :: sensible heat flux (2-dim)
61			C EVAP :: evaporation [g/(m^2 s)] (2-dim)
62			C SLRU :: sfc lw radiation (upward flux) (2-dim)
63	jmc	1.4	C Shf0 :: sensible heat flux over freezing surf.
64			C dShf :: sensible heat flux derivative relative to surf. temp
65	jmc	1.1	C Evp0 :: evaporation computed over freezing surface (Ts=0.oC)
66			C dEvp :: evaporation derivative relative to surf. temp
67			C Slr0 :: upward long wave radiation over freezing surf.
68			C dSlr :: upward long wave rad. derivative relative to surf. temp
69	jmc	1.2	C sFlx :: net heat flux (+=down) except SW, function of surf. temp Ts:
70	jmc	1.1	C 0: Flux(Ts=0.oC) ; 1: Flux(Ts^n) ; 2: d.Flux/d.Ts(Ts^n)
71			C TSFC :: surface temperature (clim.) (2-dim)
72			C TSKIN :: skin surface temperature (2-dim)
73			C-- Input:
74			C bi,bj :: tile index
75			C myThid :: Thread number for this instance of the routine
76			C--
77			_RL PSA(NGP), FMASK(NGP), EMISloc
78			_RL Tsurf(NGP), dTskin(NGP)
79			_RL SSR(NGP), SLRD(NGP)
80	jmc	1.4	_RL T1(NGP), T0(NGP), Q0(NGP), DENVV(NGP)
81	jmc	1.1
82			_RL SHF(NGP), EVAP(NGP), SLRU(NGP)
83	jmc	1.4	_RL Shf0(NGP), dShf(NGP), Evp0(NGP), dEvp(NGP)
84			_RL Slr0(NGP), dSlr(NGP), sFlx(NGP,0:2)
85	jmc	1.1	_RL TSFC(NGP), TSKIN(NGP)
86
87			INTEGER bi,bj,myThid
88			CEOP
89
90			#ifdef ALLOW_AIM
91
92			C-- Local variables:
93	jmc	1.4	C CDENVV :: surf. heat flux (sens.,lat.) coeff including stability effect
94			_RL CDENVV(NGP), RDTH, FSSICE
95			_RL Fstb0, dTstb, dFstb
96	jmc	1.1	_RL QSAT0(NGP,2)
97			_RL QDUMMY(1), RDUMMY(1), TS2
98			INTEGER J
99
100			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
101
102			C 1.5 Define effective skin temperature to compensate for
103			C non-linearity of heat/moisture fluxes during the daily cycle
104
105			DO J=1,NGP
106			c TSKIN(J) = Tsurf(J) + dTskin(J)
107			c TSFC(J)=273.16 _d 0 + dTskin(J)
108			TSKIN(J) = Tsurf(J)
109			TSFC(J)=273.16 _d 0
110			ENDDO
111
112
113			C-- 2. Computation of fluxes over land and sea
114
115	jmc	1.4	C 2.1 Stability correction
116	jmc	1.1
117	jmc	1.4	RDTH = FSTAB/DTHETA
118	jmc	1.1
119			DO J=1,NGP
120	jmc	1.4	FSSICE=1.+MIN(DTHETA,MAX(-DTHETA,TSKIN(J)-T1(J)))*RDTH
121			CDENVV(J)=CHSDENVV(J)FSSICE
122	jmc	1.1	ENDDO
123
124	jmc	1.4	IF ( dTstab.GT.0. _d 0 ) THEN
125			C- account for stability function derivative relative to Tsurf:
126			C note: to avoid discontinuity in the derivative (because of min,max), compute
127			C the derivative using the discrete form: F(Ts+dTstab)-F(Ts-dTstab)/2.dTstab
128			DO J=1,NGP
129			Fstb0 = 1.+MIN(DTHETA,MAX(-DTHETA,TSFC(J) -T1(J)))*RDTH
130			Shf0(J) = CHLDENVV(J)Fstb0
131			dTstb = ( DTHETA+dTstab-ABS(TSKIN(J)-T1(J)) )/dTstab
132			dFstb = RDTHMIN(1. _d 0, MAX(0. _d 0, dTstb0.5 _d 0))
133			dShf(J) = CHLDENVV(J)dFstb
134			ENDDO
135			ENDIF
136
137			C 2.2 Evaporation
138	jmc	1.1
139			CALL SHTORH (2, NGP, TSKIN, PSA, 1. _d 0, QDUMMY, dEvp,
140			& QSAT0(1,1), myThid)
141			CALL SHTORH (0, NGP, TSFC, PSA, 1. _d 0, QDUMMY, RDUMMY,
142			& QSAT0(1,2), myThid)
143
144	jmc	1.4	IF ( dTstab.GT.0. _d 0 ) THEN
145			C- account for stability function derivative relative to Tsurf:
146			DO J=1,NGP
147			EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
148			Evp0(J) = Shf0(J)*(QSAT0(J,2)-Q0(J))
149			dEvp(J) = CDENVV(J)*dEvp(J)
150			& + dShf(J)*(QSAT0(J,1)-Q0(J))
151			ENDDO
152			ELSE
153			DO J=1,NGP
154	jmc	1.1	EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
155			Evp0(J) = CDENVV(J)*(QSAT0(J,2)-Q0(J))
156			dEvp(J) = CDENVV(J)*dEvp(J)
157	jmc	1.4	ENDDO
158			ENDIF
159
160			C 2.3 Sensible heat flux
161
162			IF ( dTstab.GT.0. _d 0 ) THEN
163			C- account for stability function derivative relative to Tsurf:
164			DO J=1,NGP
165			SHF(J) = CDENVV(J)CP(TSKIN(J)-T0(J))
166			Shf0(J) = Shf0(J)CP(TSFC(J) -T0(J))
167			dShf(J) = CDENVV(J)*CP
168			& + dShf(J)CP(TSKIN(J)-T0(J))
169			ENDDO
170			ELSE
171			DO J=1,NGP
172			SHF(J) = CDENVV(J)CP(TSKIN(J)-T0(J))
173			Shf0(J) = CDENVV(J)CP(TSFC(J) -T0(J))
174			dShf(J) = CDENVV(J)*CP
175			dShf(J) = MAX( dShf(J), 0. _d 0 )
176			ENDDO
177			ENDIF
178	jmc	1.1
179			C 2.4 Emission of lw radiation from the surface
180
181			DO J=1,NGP
182			TS2 = TSFC(J)*TSFC(J)
183			Slr0(J) = SBCTS2TS2
184			TS2 = TSKIN(J)*TSKIN(J)
185			SLRU(J) = SBCTS2TS2
186			dSlr(J) = 4. _d 0 SBCTS2*TSKIN(J)
187			ENDDO
188
189			C-- Compute net surface heat flux and its derivative ./. surf. temp.
190			DO J=1,NGP
191	jmc	1.4	sFlx(J,0)= ( SLRD(J) - EMISloc*Slr0(J) )
192			& - ( Shf0(J) + ALHC*Evp0(J) )
193			sFlx(J,1)= ( SLRD(J) - EMISloc*SLRU(J) )
194			& - ( SHF(J) + ALHC*EVAP(J) )
195			sFlx(J,2)= -EMISloc*dSlr(J)
196			& - ( dShf(J) + ALHC*dEvp(J) )
197	jmc	1.1	ENDDO
198			IF ( aim_energPrecip ) THEN
199	jmc	1.2	C- Evap of snow/ice: substract Latent Heat of freezing from heatFlux
200	jmc	1.1	DO J=1,NGP
201			sFlx(J,0) = sFlx(J,0) - ALHF*Evp0(J)
202			sFlx(J,1) = sFlx(J,1) - ALHF*EVAP(J)
203			sFlx(J,2) = sFlx(J,2) - ALHF*dEvp(J)
204			ENDDO
205			ENDIF
206
207			C-- 3. Adjustment of skin temperature and fluxes over land
208			C-- based on energy balance (to be implemented)
209			C <= done separately for each surface type (land,ocean,sea-ice)
210
211			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
212			#endif /* ALLOW_AIM */
213
214			RETURN
215			END