C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/aim_v23/phy_suflux_sice.F,v 1.6 2006/03/13 03:58:32 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" #include "PARAMS.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 C ALHevp :: Latent Heat of evaporation _RL CDENVV(NGP), RDTH, FSSICE _RL ALHevp, Fstb0, dTstb, dFstb _RL QSAT0(NGP,2) _RL QDUMMY(1), RDUMMY(1), TS2 INTEGER J C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| ALHevp = ALHC C Evap of snow/ice: account for Latent Heat of freezing : IF ( aim_energPrecip .OR. useThSIce ) ALHevp = ALHC + ALHF 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) = CHS*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) = CHS*DENVV(J)*dFstb ENDDO C- deBug part: c J = 6 + (17-1)*sNx c IF ( bi.EQ.3 .AND. J.LE.NGP ) c & WRITE(6,1020)'SUFLUX_SICE: Stab=',Shf0(J),CDENVV(J),dShf(J) 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)) dShf(J) = MAX( dShf(J), 0. _d 0 ) C-- do not allow negative derivative vs Ts of Sensible+Latent H.flux: C a) quiet unrealistic ; C b) garantee positive deriv. of total H.flux (needed for implicit solver) dEvp(J) = MAX( dEvp(J), -dShf(J)/ALHevp ) 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 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) + ALHevp*Evp0(J) ) sFlx(J,1)= ( SLRD(J) - EMISloc*SLRU(J) ) & - ( SHF(J) + ALHevp*EVAP(J) ) sFlx(J,2)= -EMISloc*dSlr(J) & - ( dShf(J) + ALHevp*dEvp(J) ) ENDDO C- deBug part: ----------------- c1010 FORMAT(A,I3,2F10.3,F10.4) c1020 FORMAT(A,1P4E11.3) c J = 6 + (17-1)*sNx c IF ( bi.EQ.3 .AND. J.LE.NGP ) THEN c WRITE(6,1010) 'SUFLUX_SICE: 1,sFlx=', 1, c & sFlx(J,0),sFlx(J,1),sFlx(J,2) c WRITE(6,1010) 'SUFLUX_SICE: 0,Evap=', 0,Evp0(J),EVAP(J),dEvp(J) c WRITE(6,1010) 'SUFLUX_SICE: -,LWup=',-1,Slr0(J),SLRU(J),dSlr(J) c WRITE(6,1010) 'SUFLUX_SICE: -, SHF=',-1,Shf0(J),SHF(J), dShf(J) c WRITE(6,1010) 'SUFLUX_SICE: -, LAT=',-1, c & ALHevp*Evp0(J),ALHevp*EVAP(J),ALHevp*dEvp(J) c 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