Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:38:12 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
88ac963ab8 Jean*0001 #include "BULK_FORCE_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: BULKF_FORMULA_AIM
                0005 C     !INTERFACE:
                0006       SUBROUTINE BULKF_FORMULA_AIM(
                0007      I                   Tsurf, SLRD,
                0008      I                   T1, T0, Q0, Vsurf,
                0009      O                   SHF, EVAP, SLRU,
                0010      O                   dEvp, sFlx,
                0011      I                   iceornot, myThid )
                0012 
                0013 C     !DESCRIPTION: \bv
                0014 C     *==========================================================*
                0015 C     | S/R BULKF_FORMULA_AIM
                0016 C     | o compute surface flux over ocean and sea-ice,
                0017 C     |   using AIM surface flux formulation
                0018 C     *==========================================================*
                0019 C     *==========================================================*
                0020 C     \ev
                0021 
                0022 C     !USES:
                0023       IMPLICIT NONE
                0024 
                0025 C     Resolution parameters
                0026 
                0027 #include "EEPARAMS.h"
                0028 #include "SIZE.h"
                0029 #include "PARAMS.h"
                0030 #include "BULKF_PARAMS.h"
                0031 
                0032 C     !INPUT/OUTPUT PARAMETERS:
                0033 C     == Routine Arguments ==
                0034 C--   Input:
                0035 C    FMASK  :: fractional land-sea mask        (2-dim)
                0036 C    Tsurf  :: surface temperature        (2-dim)
                0037 C    SSR    :: sfc sw radiation (net flux)     (2-dim)
                0038 C    SLRD   :: sfc lw radiation (downward flux)(2-dim)
                0039 C    T1     :: near-surface air temperature (from Pot.temp)
                0040 C    T0     :: near-surface air temperature    (2-dim)
                0041 C    Q0     :: near-surface sp. humidity [g/kg](2-dim)
                0042 C    Vsurf  :: surface wind speed        [m/s] (2-dim,input)
                0043 C--   Output:
                0044 C    SHF    :: sensible heat flux              (2-dim)
                0045 C    EVAP   :: evaporation [g/(m^2 s)]         (2-dim)
                0046 C    SLRU   :: sfc lw radiation (upward flux)  (2-dim)
                0047 C    Shf0   :: sensible heat flux over freezing surf.
                0048 C    dShf   :: sensible heat flux derivative relative to surf. temp
                0049 C    Evp0   :: evaporation computed over freezing surface (Ts=0.oC)
                0050 C    dEvp   :: evaporation derivative relative to surf. temp
                0051 C    Slr0   :: upward long wave radiation over freezing surf.
                0052 C    dSlr   :: upward long wave rad. derivative relative to surf. temp
                0053 C    sFlx   :: net heat flux (+=down) except SW, function of surf. temp Ts:
                0054 C              0: Flux(Ts=0.oC) ; 1: Flux(Ts^n) ; 2: d.Flux/d.Ts(Ts^n)
                0055 C    TSFC   :: surface temperature (clim.)     (2-dim)
                0056 C    TSKIN  :: skin surface temperature        (2-dim)
                0057 C--   Input:
                0058 C    iceornot :: 0=open water, 1=ice cover
                0059 C    myThid :: Thread number for this instance of the routine
                0060 C--
                0061       INTEGER NGP
                0062       PARAMETER ( NGP = 1 )
                0063 c     _RL  PSA(NGP), FMASK(NGP), EMISloc
                0064       _RL  Tsurf(NGP)
                0065 c     _RL  SSR(NGP)
                0066       _RL  SLRD(NGP)
                0067       _RL  T1(NGP), T0(NGP), Q0(NGP), Vsurf(NGP)
                0068 
                0069       _RL  SHF(NGP), EVAP(NGP), SLRU(NGP)
                0070       _RL  dEvp(NGP), sFlx(NGP,0:2)
                0071 c     _RL  Shf0(NGP), dShf(NGP), Evp0(NGP), dEvp(NGP)
                0072 c     _RL  Slr0(NGP), dSlr(NGP), sFlx(NGP,0:2)
                0073 c     _RL  TSFC(NGP), TSKIN(NGP)
                0074 
                0075       INTEGER iceornot
                0076       INTEGER myThid
                0077 CEOP
                0078 
                0079 #ifdef ALLOW_FORMULA_AIM
                0080 
                0081 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0082 C      FWIND0 = ratio of near-sfc wind to lowest-level wind
                0083 C      CHS    = heat exchange coefficient over sea
                0084 C      VGUST  = wind speed for sub-grid-scale gusts
                0085 C      DTHETA = Potential temp. gradient for stability correction
                0086 C      dTstab = potential temp. increment for stability function derivative
                0087 C      FSTAB  = Amplitude of stability correction (fraction)
                0088 C       P0    = reference pressure                 [Pa=N/m2]
                0089 C       GG    = gravity accel.                     [m/s2]
                0090 C       RD    = gas constant for dry air           [J/kg/K]
                0091 C       CP    = specific heat at constant pressure [J/kg/K]
                0092 C       ALHC  = latent heat of condensation        [J/g]
                0093 C       ALHF  = latent heat of freezing            [J/g]
                0094 C       SBC   = Stefan-Boltzmann constant
                0095 C      EMISloc :: longwave surface emissivity
                0096 c     _RL FWIND0, CHS, VGUST, DTHETA, dTstab, FSTAB
c8e7a08f49 Jean*0097       _RL P0, ALHC, ALHF, RD, CP, SBC, EMISloc
88ac963ab8 Jean*0098       EQUIVALENCE ( ocean_emissivity , EMISloc )
                0099       EQUIVALENCE ( Lvap   , ALHC )
                0100       EQUIVALENCE ( Lfresh , ALHF )
                0101       EQUIVALENCE ( Rgas   , RD )
                0102       EQUIVALENCE ( cpair  , CP )
                0103       EQUIVALENCE ( stefan , SBC )
                0104 
                0105 C-- Local variables:
                0106 C    PSA     :: norm. surface pressure [p/p0]   (2-dim)
                0107 C    DENVV   :: surface flux (sens,lat.) coeff. (=Rho*|V|) [kg/m2/s]
                0108 C    CDENVV  :: surf. heat flux (sens.,lat.) coeff including stability effect
                0109 C    ALHevp  :: Latent Heat of evaporation
                0110       _RL PSA(NGP)
c8e7a08f49 Jean*0111       _RL DENVV(NGP), PRD
88ac963ab8 Jean*0112       _RL  Shf0(NGP), dShf(NGP), Evp0(NGP)
                0113       _RL  Slr0(NGP), dSlr(NGP)
                0114       _RL  TSFC(NGP), TSKIN(NGP)
                0115       _RL CDENVV(NGP), RDTH, FSSICE
                0116       _RL ALHevp, Fstb0, dTstb, dFstb
                0117       _RL QSAT0(NGP,2)
                0118       _RL QDUMMY(1), RDUMMY(1), TS2
                0119       INTEGER J
                0120 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0121 
                0122       PSA(1) = 1. _d 0
c8e7a08f49 Jean*0123       P0 = 1. _d +5
88ac963ab8 Jean*0124 C---
                0125 
                0126       ALHevp = ALHC
                0127 C     Evap of snow/ice: account for Latent Heat of freezing :
c8e7a08f49 Jean*0128 c     IF ( aim_energPrecip .OR. useThSIce ) ALHevp = ALHC + ALHF
                0129       IF ( iceornot.GE.1 ) ALHevp = ALHC + ALHF
88ac963ab8 Jean*0130 
                0131 C     1.4 Density * wind speed (including gustiness factor)
                0132 
                0133       PRD = P0/RD
c8e7a08f49 Jean*0134 c     VG2 = VGUST*VGUST
                0135 c     factWind2 = FWIND0*FWIND0
88ac963ab8 Jean*0136 
                0137       DO J=1,NGP
                0138 c       SPEED0(J)=SQRT(factWind2*Vsurf2(J)+VG2)
                0139 c       DENVV(J)=(PRD*PSA(J)/T0(J))*SPEED0(J)
c8e7a08f49 Jean*0140 C-- assuming input file "WspeedFile" contains the time-average "SPEED0"
                0141 C     from AIM output (aimPhytave: fields # 15 ; aimDiag: WINDS ) :
88ac963ab8 Jean*0142         DENVV(J)=(PRD*PSA(J)/T0(J))*Vsurf(J)
                0143       ENDDO
                0144 
                0145 C     1.5 Define effective skin temperature to compensate for
                0146 C         non-linearity of heat/moisture fluxes during the daily cycle
                0147 
                0148       DO J=1,NGP
                0149         TSKIN(J) = Tsurf(J) + celsius2K
                0150         TSFC(J)=273.16 _d 0
                0151       ENDDO
                0152 
                0153 C--   2. Computation of fluxes over land and sea
                0154 
                0155 C     2.1 Stability correction
                0156 
                0157       RDTH = FSTAB/DTHETA
                0158 
                0159       DO J=1,NGP
                0160         FSSICE=1.+MIN(DTHETA,MAX(-DTHETA,TSKIN(J)-T1(J)))*RDTH
                0161         CDENVV(J)=CHS*DENVV(J)*FSSICE
                0162       ENDDO
                0163 
                0164       IF ( dTstab.GT.0. _d 0 ) THEN
                0165 C-    account for stability function derivative relative to Tsurf:
                0166 C note: to avoid discontinuity in the derivative (because of min,max), compute
                0167 C   the derivative using the discrete form: F(Ts+dTstab)-F(Ts-dTstab)/2.dTstab
                0168        DO J=1,NGP
                0169         Fstb0 = 1.+MIN(DTHETA,MAX(-DTHETA,TSFC(J) -T1(J)))*RDTH
                0170         Shf0(J) = CHS*DENVV(J)*Fstb0
                0171         dTstb = ( DTHETA+dTstab-ABS(TSKIN(J)-T1(J)) )/dTstab
                0172         dFstb = RDTH*MIN(1. _d 0, MAX(0. _d 0, dTstb*0.5 _d 0))
                0173         dShf(J) = CHS*DENVV(J)*dFstb
                0174        ENDDO
                0175       ENDIF
                0176 
                0177 C     2.2 Evaporation
                0178 
                0179       CALL BULKF_SH2RH_AIM( 2, NGP, TSKIN, PSA, 1. _d 0,
                0180      &                      QDUMMY, dEvp,  QSAT0(1,1), myThid )
                0181       CALL BULKF_SH2RH_AIM( 0, NGP, TSFC,  PSA, 1. _d 0,
                0182      &                      QDUMMY, RDUMMY,QSAT0(1,2), myThid )
                0183 
                0184       IF ( dTstab.GT.0. _d 0 ) THEN
                0185 C-    account for stability function derivative relative to Tsurf:
                0186        DO J=1,NGP
                0187         EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
                0188         Evp0(J) =   Shf0(J)*(QSAT0(J,2)-Q0(J))
                0189         dEvp(J) = CDENVV(J)*dEvp(J)
                0190      &            + dShf(J)*(QSAT0(J,1)-Q0(J))
                0191        ENDDO
                0192       ELSE
                0193        DO J=1,NGP
                0194         EVAP(J) = CDENVV(J)*(QSAT0(J,1)-Q0(J))
                0195         Evp0(J) = CDENVV(J)*(QSAT0(J,2)-Q0(J))
                0196         dEvp(J) = CDENVV(J)*dEvp(J)
                0197        ENDDO
                0198       ENDIF
                0199 
                0200 C     2.3 Sensible heat flux
                0201 
                0202       IF ( dTstab.GT.0. _d 0 ) THEN
                0203 C-    account for stability function derivative relative to Tsurf:
                0204        DO J=1,NGP
                0205         SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
                0206         Shf0(J) =   Shf0(J)*CP*(TSFC(J) -T0(J))
                0207         dShf(J) = CDENVV(J)*CP
                0208      &            + dShf(J)*CP*(TSKIN(J)-T0(J))
                0209         dShf(J) = MAX( dShf(J), 0. _d 0 )
                0210 C--   do not allow negative derivative vs Ts of Sensible+Latent H.flux:
                0211 C     a) quiet unrealistic ;
                0212 C     b) garantee positive deriv. of total H.flux (needed for implicit solver)
                0213         dEvp(J) = MAX( dEvp(J), -dShf(J)/ALHevp )
                0214        ENDDO
                0215       ELSE
                0216        DO J=1,NGP
                0217         SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))
                0218         Shf0(J) = CDENVV(J)*CP*(TSFC(J) -T0(J))
                0219         dShf(J) = CDENVV(J)*CP
                0220        ENDDO
                0221       ENDIF
                0222 
                0223 C     2.4 Emission of lw radiation from the surface
                0224 
                0225       DO J=1,NGP
                0226         TS2     = TSFC(J)*TSFC(J)
                0227         Slr0(J) = SBC*TS2*TS2
                0228         TS2     = TSKIN(J)*TSKIN(J)
                0229         SLRU(J) = SBC*TS2*TS2
                0230         dSlr(J)  = 4. _d 0 *SBC*TS2*TSKIN(J)
                0231       ENDDO
                0232 
                0233 C--   Compute net surface heat flux and its derivative ./. surf. temp.
                0234       DO J=1,NGP
                0235         sFlx(J,0)= ( SLRD(J) - EMISloc*Slr0(J) )
                0236      &           - ( Shf0(J) + ALHevp*Evp0(J) )
                0237         sFlx(J,1)= ( SLRD(J) - EMISloc*SLRU(J) )
                0238      &           - ( SHF(J)  + ALHevp*EVAP(J) )
                0239         sFlx(J,2)=            -EMISloc*dSlr(J)
                0240      &           - ( dShf(J) + ALHevp*dEvp(J) )
                0241       ENDDO
                0242 
                0243 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0244 #endif /* ALLOW_FORMULA_AIM */
                0245 
                0246       RETURN
                0247       END