pkg/aim_v23/phy_suflux_sice.F

b3097ed02d Jean*0001 #include "AIM_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SUFLUX_SICE
                0005 C     !INTERFACE:
                0006       SUBROUTINE SUFLUX_SICE(
                0007      I                   PSA, FMASK, EMISloc,
                0008      I                   Tsurf, dTskin, SSR, SLRD,
e749d70ece Jean*0009      I                   T1, T0, Q0, DENVV,
b3097ed02d Jean*0010      O                   SHF, EVAP, SLRU,
e749d70ece Jean*0011      O                   Shf0, dShf, Evp0, dEvp, Slr0, dSlr, sFlx,
b3097ed02d Jean*0012      O                   TSFC, TSKIN,
                0013      I                   bi,bj,myThid)
                0014
                0015 C     !DESCRIPTION: \bv
                0016 C     *==========================================================*
                0017 C     | S/R SUFLUX_SICE
                0018 C     | o compute surface flux over sea-ice
                0019 C     *==========================================================*
                0020 C     | o contains part of original S/R SUFLUX (Speedy code)
                0021 C     *==========================================================*
                0022 C     \ev
                0023
                0024 C     !USES:
                0025       IMPLICIT NONE
                0026
                0027 C     Resolution parameters
                0028
                0029 C-- size for MITgcm & Physics package :
                0030 #include "AIM_SIZE.h"
                0031 #include "EEPARAMS.h"
a232c4b875 Jean*0032 #include "PARAMS.h"
b3097ed02d Jean*0033
                0034 C-- Physics package
                0035 #include "AIM_PARAMS.h"
                0036
                0037 C     Physical constants + functions of sigma and latitude
                0038 #include "com_physcon.h"
                0039
                0040 C     Surface flux constants
                0041 #include "com_sflcon.h"
                0042
                0043 C     !INPUT/OUTPUT PARAMETERS:
                0044 C     == Routine Arguments ==
                0045 C--   Input:
                0046 C    PSA    :: norm. surface pressure [p/p0]   (2-dim)
                0047 C    FMASK  :: fractional land-sea mask        (2-dim)
                0048 C    EMISloc:: longwave surface emissivity
                0049 C    Tsurf  :: surface temperature        (2-dim)
                0050 C    dTskin :: temp. correction for daily-cycle heating [K]
                0051 C    SSR    :: sfc sw radiation (net flux)     (2-dim)
                0052 C    SLRD   :: sfc lw radiation (downward flux)(2-dim)
e749d70ece Jean*0053 C    T1     :: near-surface air temperature (from Pot.temp)
b3097ed02d Jean*0054 C    T0     :: near-surface air temperature    (2-dim)
                0055 C    Q0     :: near-surface sp. humidity [g/kg](2-dim)
e749d70ece Jean*0056 C    DENVV  :: surface flux (sens,lat.) coeff. (=Rho*|V|) [kg/m2/s]
b3097ed02d Jean*0057 C--   Output:
                0058 C    SHF    :: sensible heat flux              (2-dim)
                0059 C    EVAP   :: evaporation [g/(m^2 s)]         (2-dim)
                0060 C    SLRU   :: sfc lw radiation (upward flux)  (2-dim)
e749d70ece Jean*0061 C    Shf0   :: sensible heat flux over freezing surf.
                0062 C    dShf   :: sensible heat flux derivative relative to surf. temp
b3097ed02d Jean*0063 C    Evp0   :: evaporation computed over freezing surface (Ts=0.oC)
                0064 C    dEvp   :: evaporation derivative relative to surf. temp
                0065 C    Slr0   :: upward long wave radiation over freezing surf.
                0066 C    dSlr   :: upward long wave rad. derivative relative to surf. temp
cdcb187d4c Jean*0067 C    sFlx   :: net heat flux (+=down) except SW, function of surf. temp Ts:
b3097ed02d Jean*0068 C              0: Flux(Ts=0.oC) ; 1: Flux(Ts^n) ; 2: d.Flux/d.Ts(Ts^n)
                0069 C    TSFC   :: surface temperature (clim.)     (2-dim)
                0070 C    TSKIN  :: skin surface temperature        (2-dim)
                0071 C--   Input:
                0072 C    bi,bj  :: tile index
                0073 C    myThid :: Thread number for this instance of the routine
                0074 C--
                0075       _RL  PSA(NGP), FMASK(NGP), EMISloc
                0076       _RL  Tsurf(NGP), dTskin(NGP)
                0077       _RL  SSR(NGP), SLRD(NGP)
e749d70ece Jean*0078       _RL  T1(NGP), T0(NGP), Q0(NGP), DENVV(NGP)
b3097ed02d Jean*0079
                0080       _RL  SHF(NGP), EVAP(NGP), SLRU(NGP)
e749d70ece Jean*0081       _RL  Shf0(NGP), dShf(NGP), Evp0(NGP), dEvp(NGP)
                0082       _RL  Slr0(NGP), dSlr(NGP), sFlx(NGP,0:2)
b3097ed02d Jean*0083       _RL  TSFC(NGP), TSKIN(NGP)
                0084
                0085       INTEGER bi,bj,myThid
                0086 CEOP
                0087
                0088 #ifdef ALLOW_AIM
                0089
                0090 C-- Local variables:
e749d70ece Jean*0091 C    CDENVV :: surf. heat flux (sens.,lat.) coeff including stability effect
a232c4b875 Jean*0092 C    ALHevp :: Latent Heat of evaporation
e749d70ece Jean*0093       _RL CDENVV(NGP), RDTH, FSSICE
097d65e6b5 Jean*0094       _RL ALHevp, Fstb0, dTstb, dFstb
b3097ed02d Jean*0095       _RL QSAT0(NGP,2)
                0096       _RL QDUMMY(1), RDUMMY(1), TS2
                0097       INTEGER J
                0098
                0099 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0100
097d65e6b5 Jean*0101       ALHevp = ALHC
                0102 C     Evap of snow/ice: account for Latent Heat of freezing :
a232c4b875 Jean*0103       IF ( aim_energPrecip .OR. useThSIce ) ALHevp = ALHC + ALHF
097d65e6b5 Jean*0104
b3097ed02d Jean*0105 C     1.5 Define effective skin temperature to compensate for
                0106 C         non-linearity of heat/moisture fluxes during the daily cycle
                0107
                0108       DO J=1,NGP
                0109 c       TSKIN(J) = Tsurf(J) + dTskin(J)
                0110 c       TSFC(J)=273.16 _d 0 + dTskin(J)
                0111         TSKIN(J) = Tsurf(J)
                0112         TSFC(J)=273.16 _d 0
                0113       ENDDO
                0114
                0115 C--   2. Computation of fluxes over land and sea
                0116
e749d70ece Jean*0117 C     2.1 Stability correction
b3097ed02d Jean*0118
e749d70ece Jean*0119       RDTH = FSTAB/DTHETA
b3097ed02d Jean*0120
                0121       DO J=1,NGP
e749d70ece Jean*0122         FSSICE=1.+MIN(DTHETA,MAX(-DTHETA,TSKIN(J)-T1(J)))*RDTH
                0123         CDENVV(J)=CHS*DENVV(J)*FSSICE
b3097ed02d Jean*0124       ENDDO
                0125
e749d70ece Jean*0126       IF ( dTstab.GT.0. _d 0 ) THEN
                0127 C-    account for stability function derivative relative to Tsurf:
                0128 C note: to avoid discontinuity in the derivative (because of min,max), compute
                0129 C   the derivative using the discrete form: F(Ts+dTstab)-F(Ts-dTstab)/2.dTstab
                0130        DO J=1,NGP
                0131         Fstb0 = 1.+MIN(DTHETA,MAX(-DTHETA,TSFC(J) -T1(J)))*RDTH
097d65e6b5 Jean*0132         Shf0(J) = CHS*DENVV(J)*Fstb0
e749d70ece Jean*0133         dTstb = ( DTHETA+dTstab-ABS(TSKIN(J)-T1(J)) )/dTstab
                0134         dFstb = RDTH*MIN(1. _d 0, MAX(0. _d 0, dTstb*0.5 _d 0))
097d65e6b5 Jean*0135         dShf(J) = CHS*DENVV(J)*dFstb
e749d70ece Jean*0136        ENDDO
097d65e6b5 Jean*0137 C-    deBug part:
                0138 c      J = 6 + (17-1)*sNx
a232c4b875 Jean*0139 c      IF ( bi.EQ.3 .AND. J.LE.NGP )
097d65e6b5 Jean*0140 c    &  WRITE(6,1020)'SUFLUX_SICE: Stab=',Shf0(J),CDENVV(J),dShf(J)
e749d70ece Jean*0141       ENDIF
                0142
                0143 C     2.2 Evaporation
b3097ed02d Jean*0144
                0145       CALL SHTORH (2, NGP, TSKIN, PSA, 1. _d 0, QDUMMY, dEvp,
                0146      &                QSAT0(1,1), myThid)
                0147       CALL SHTORH (0, NGP, TSFC, PSA, 1. _d 0, QDUMMY, RDUMMY,
                0148      &                QSAT0(1,2), myThid)
                0149
e749d70ece Jean*0150       IF ( dTstab.GT.0. _d 0 ) THEN
                0151 C-    account for stability function derivative relative to Tsurf:
                0152        DO J=1,NGP
                0153         EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
                0154         Evp0(J) =   Shf0(J)*(QSAT0(J,2)-Q0(J))
                0155         dEvp(J) = CDENVV(J)*dEvp(J)
                0156      &            + dShf(J)*(QSAT0(J,1)-Q0(J))
                0157        ENDDO
                0158       ELSE
                0159        DO J=1,NGP
b3097ed02d Jean*0160         EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
                0161         Evp0(J) = CDENVV(J)*(QSAT0(J,2)-Q0(J))
                0162         dEvp(J) = CDENVV(J)*dEvp(J)
e749d70ece Jean*0163        ENDDO
                0164       ENDIF
                0165
                0166 C     2.3 Sensible heat flux
                0167
                0168       IF ( dTstab.GT.0. _d 0 ) THEN
                0169 C-    account for stability function derivative relative to Tsurf:
                0170        DO J=1,NGP
                0171         SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
                0172         Shf0(J) =   Shf0(J)*CP*(TSFC(J) -T0(J))
                0173         dShf(J) = CDENVV(J)*CP
                0174      &            + dShf(J)*CP*(TSKIN(J)-T0(J))
097d65e6b5 Jean*0175         dShf(J) = MAX( dShf(J), 0. _d 0 )
                0176 C--   do not allow negative derivative vs Ts of Sensible+Latent H.flux:
a232c4b875 Jean*0177 C     a) quiet unrealistic ;
097d65e6b5 Jean*0178 C     b) garantee positive deriv. of total H.flux (needed for implicit solver)
                0179         dEvp(J) = MAX( dEvp(J), -dShf(J)/ALHevp )
e749d70ece Jean*0180        ENDDO
                0181       ELSE
                0182        DO J=1,NGP
                0183         SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
                0184         Shf0(J) = CDENVV(J)*CP*(TSFC(J) -T0(J))
                0185         dShf(J) = CDENVV(J)*CP
                0186        ENDDO
                0187       ENDIF
b3097ed02d Jean*0188
                0189 C     2.4 Emission of lw radiation from the surface
                0190
                0191       DO J=1,NGP
                0192         TS2     = TSFC(J)*TSFC(J)
                0193         Slr0(J) = SBC*TS2*TS2
                0194         TS2     = TSKIN(J)*TSKIN(J)
                0195         SLRU(J) = SBC*TS2*TS2
                0196         dSlr(J)  = 4. _d 0 *SBC*TS2*TSKIN(J)
                0197       ENDDO
                0198
                0199 C--   Compute net surface heat flux and its derivative ./. surf. temp.
                0200       DO J=1,NGP
e749d70ece Jean*0201         sFlx(J,0)= ( SLRD(J) - EMISloc*Slr0(J) )
097d65e6b5 Jean*0202      &           - ( Shf0(J) + ALHevp*Evp0(J) )
e749d70ece Jean*0203         sFlx(J,1)= ( SLRD(J) - EMISloc*SLRU(J) )
097d65e6b5 Jean*0204      &           - ( SHF(J)  + ALHevp*EVAP(J) )
e749d70ece Jean*0205         sFlx(J,2)=            -EMISloc*dSlr(J)
097d65e6b5 Jean*0206      &           - ( dShf(J) + ALHevp*dEvp(J) )
b3097ed02d Jean*0207       ENDDO
097d65e6b5 Jean*0208
                0209 C-    deBug part:   -----------------
                0210 c1010 FORMAT(A,I3,2F10.3,F10.4)
                0211 c1020 FORMAT(A,1P4E11.3)
                0212 c     J = 6 + (17-1)*sNx
                0213 c     IF ( bi.EQ.3 .AND. J.LE.NGP ) THEN
                0214 c       WRITE(6,1010) 'SUFLUX_SICE: 1,sFlx=', 1,
                0215 c    &                                    sFlx(J,0),sFlx(J,1),sFlx(J,2)
                0216 c       WRITE(6,1010) 'SUFLUX_SICE: 0,Evap=', 0,Evp0(J),EVAP(J),dEvp(J)
                0217 c       WRITE(6,1010) 'SUFLUX_SICE: -,LWup=',-1,Slr0(J),SLRU(J),dSlr(J)
                0218 c       WRITE(6,1010) 'SUFLUX_SICE: -, SHF=',-1,Shf0(J),SHF(J), dShf(J)
                0219 c       WRITE(6,1010) 'SUFLUX_SICE: -, LAT=',-1,
                0220 c    &                     ALHevp*Evp0(J),ALHevp*EVAP(J),ALHevp*dEvp(J)
                0221 c     ENDIF
b3097ed02d Jean*0222
                0223 C--   3. Adjustment of skin temperature and fluxes over land
                0224 C--      based on energy balance (to be implemented)
                0225 C        <= done separately for each surface type (land,ocean,sea-ice)
                0226
                0227 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0228 #endif /* ALLOW_AIM */
                0229
                0230       RETURN
                0231       END