Changeset 11772 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
- Timestamp:
- 2019-10-23T16:04:12+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
r11637 r11772 38 38 39 39 !! Cool-skin related parameters: 40 ! ... 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 41 & dT_cs !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 42 REAL(wp), PARAMETER :: zcon0 = -16._wp * grav * rho0_w * rCp0_w * rnu0_w*rnu0_w*rnu0_w / ( rk0_w*rk0_w ) ! "-" because ocean convention: Qabs > 0 => gain of heat for ocean! 43 !! => see eq.(14) in Fairall et al. 1996 (eq.(6) of Zeng aand Beljaars is WRONG! (typo?) 41 44 42 45 !! Warm-layer related parameters: 43 REAL(wp), PARAMETER :: rd0 = 3. !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) 44 REAL(wp), PARAMETER :: zRhoCp_w = rho0_w*rCp0_w 45 REAL(wp), PARAMETER :: rNu0 = 1.0 !: be closer to COARE3p6 ???!LOLO 46 !REAL(wp), PARAMETER :: rNu0 = 0.5 !: Nu (exponent of temperature profile) Eq.11 47 ! !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 48 ! !: 0.3 to respect a warming of +3 K in calm 49 ! !: condition for the insolation peak of +1000W/m^2 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 47 & dT_wl, & !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 48 & Hz_wl !: depth of warm-layer [m] 49 ! 50 REAL(wp), PARAMETER, PUBLIC :: rd0 = 3. !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) 51 REAL(wp), PARAMETER :: zRhoCp_w = rho0_w*rCp0_w 52 ! 53 REAL(wp), PARAMETER :: rNuwl0 = 0.5 !: Nu (exponent of temperature profile) Eq.11 54 ! !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 55 ! !: 0.3 to respect a warming of +3 K in calm 56 ! !: condition for the insolation peak of +1000W/m^2 50 57 !!---------------------------------------------------------------------- 51 58 CONTAINS 52 59 53 60 54 SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST, pdT ) 55 !!--------------------------------------------------------------------- 56 !! 57 !! Cool-Skin scheme according to Fairall et al. 1996 58 !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 59 !! " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 60 !! as parameterized in IFS Cy45r1 / E.C.M.W.F. 61 !! ------------------------------------------------------------------ 61 SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST ) 62 !!--------------------------------------------------------------------- 63 !! 64 !! Cool-skin parameterization, based on Fairall et al., 1996: 65 !! 66 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 67 !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 68 !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 69 !! doi:10.1029/95JC03190. 70 !! 71 !!------------------------------------------------------------------ 62 72 !! 63 73 !! ** INPUT: … … 66 76 !! *pustar* friction velocity u* [m/s] 67 77 !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] 68 !!69 !! ** INPUT/OUTPUT:70 !! *pdT* : as input => previous estimate of dT cool-skin71 !! as output => new estimate of dT cool-skin72 !!73 78 !!------------------------------------------------------------------ 74 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 75 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 76 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) 77 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 78 !! 79 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT ! dT due to cool-skin effect 80 !!--------------------------------------------------------------------- 81 INTEGER :: ji, jj ! dummy loop indices 82 REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, & 83 & zusw, zusw2 84 !!--------------------------------------------------------------------- 85 79 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 80 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 81 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) 82 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 83 !!--------------------------------------------------------------------- 84 INTEGER :: ji, jj, jc 85 REAL(wp) :: zQabs, zdelta, zfr 86 !!--------------------------------------------------------------------- 86 87 DO jj = 1, jpj 87 88 DO ji = 1, jpi 88 89 89 zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 90 91 zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere 92 93 zusw = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw) ! u* in the water 94 zusw2 = zusw*zusw 95 96 zlamb = 6._wp*( 1._wp + (zQnsol*zalpha_w*rcst_cs/(zusw2*zusw2 ))**0.75 )**(-1./3.) ! w.r.t COARE 3p6 => seems to ommit absorbed zfr*Qsw (Qnet i.o. Qnsol) and effect of evap... 97 ! ! so zlamb not implicit in terms of zdelta (zfr(delta)), so no need to have last guess of delta as in COARE 3.6... 98 zdelta = zlamb*rnu0_w/zusw 99 100 zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq. 8.131 / IFS cy40r1, doc, Part IV, 101 zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface 102 103 !! Update! 104 pdT(ji,jj) = MIN( - zQnet*zdelta/rk0_w , 0._wp ) ! temperature increment 90 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$ 91 ! ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$ 92 !zQabs = pQnsol(ji,jj) 93 94 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 95 96 DO jc = 1, 4 ! because implicit in terms of zdelta... 97 zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 98 ! => (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 99 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface 100 !zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 101 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 102 END DO 103 104 dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp ) ! temperature increment 105 105 106 106 END DO … … 111 111 112 112 113 SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pdT)113 SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pustk ) 114 114 !!--------------------------------------------------------------------- 115 115 !! 116 116 !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 117 117 !! " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 118 !! 119 !! STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER! 118 120 !! 119 121 !! As included in IFS Cy45r1 / E.C.M.W.F. … … 125 127 !! *pustar* friction velocity u* [m/s] 126 128 !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] 127 !!128 !! ** INPUT/OUTPUT:129 !! *pdT* : as input => previous estimate of dT warm-layer130 !! as output => new estimate of dT warm-layer131 !!132 129 !!------------------------------------------------------------------ 133 130 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! surface net solar radiation into the ocean [W/m^2] => >= 0 ! … … 135 132 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity [m/s] 136 133 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST at depth gdept_1d(1) [K] 137 ! 138 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST139 ! 140 INTEGER :: ji, jj134 !! 135 REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pustk ! surface Stokes velocity [m/s] 136 ! 137 INTEGER :: ji, jj, jc 141 138 ! 142 139 REAL(wp) :: & 143 & zdz, & !: thickness of the warm-layer [m] 144 & zalpha_w, & !: thermal expansion coefficient of sea-water 145 & ZSRD, & 146 & dT_wl, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 147 & zfr,zdL,zdL2, ztmp, & 148 & ZPHI, & 149 & zus_a, zusw, zusw2, & 150 & flg, zQabs, ZL1, ZL2 151 !!--------------------------------------------------------------------- 140 & zHwl, & !: thickness of the warm-layer [m] 141 & ztcorr, & !: correction of dT w.r.t measurement depth of bulk SST (first T-point) 142 & zalpha, & !: thermal expansion coefficient of sea-water [1/K] 143 & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 144 & zfr, zeta, ztmp, & 145 & zusw, zusw2, & 146 & zLa, zfLa, & 147 & flg, zwf, zQabs, & 148 & ZA, ZB, zL1, zL2, & 149 & zcst0, zcst1, zcst2, zcst3 150 ! 151 LOGICAL :: l_pustk_known 152 !!--------------------------------------------------------------------- 153 154 l_pustk_known = .FALSE. 155 IF ( PRESENT(pustk) ) l_pustk_known = .TRUE. 152 156 153 157 DO jj = 1, jpj 154 158 DO ji = 1, jpi 155 159 156 zdz = rd0 ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !) 157 158 ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz) 160 zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !) 161 ! it is = rd0 (3m) in default Zeng & Beljaars case... 162 163 !! Previous value of dT / warm-layer, adapted to depth: 164 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$ 165 ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl 166 zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp ) 167 ! zdTwl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zHwl) 159 168 ! pdT " " and depth of bulk SST (here gdept_1d(1))! 160 !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced! 161 !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof dT_wl ! 162 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz ) ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$ 163 dT_wl = pdT(ji,jj) / ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 164 165 zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 169 !! => but of course in general the bulk SST is taken shallower than zHwl !!! So correction less pronounced! 170 !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl ! 171 172 zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 166 173 167 174 168 175 ! *** zfr = Fraction of solar radiation absorbed in warm layer (-) 169 zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*z dz) - 0.27_wp*EXP(-2.8_wp*zdz) - 0.45_wp*EXP(-0.07_wp*zdz) !: Eq. 8.157176 zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*zHwl) - 0.27_wp*EXP(-2.8_wp*zHwl) - 0.45_wp*EXP(-0.07_wp*zHwl) !: Eq. 8.157 170 177 171 178 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! tot heat absorbed in warm layer 172 179 173 zusw = MAX( pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)! u* in the water180 zusw = MAX( pustar(ji,jj), 1.E-4_wp ) * sq_radrw ! u* in the water 174 181 zusw2 = zusw*zusw 175 182 176 177 !! *** 1st rhs term in eq. 8.156 (IFS doc Cy45r1): 178 ZL1 = zQabs / ( zdz * zRhoCp_w * rNu0 ) * (rNu0 + 1._wp) 179 180 181 !! Buoyancy flux and stability parameter (zdl = -z/L) in water 182 ZSRD = zQabs/zRhoCp_w 183 ! 184 flg = 0.5_wp + SIGN(0.5_wp, ZSRD) ! ZSRD > 0. => 1. / ZSRD < 0. => 0. 185 ztmp = MAX(dT_wl,0._wp) 186 zdl = (flg+1._wp) * ( zusw2 * SQRT(ztmp/(5._wp*zdz*grav*zalpha_w/rNu0)) ) & ! (dT_wl > 0.0 .AND. ZSRD < 0.0) 187 & + flg * ZSRD ! otherwize 188 ! 189 zus_a = MAX( pustar(ji,jj), 1.E-4_wp ) 190 zdL = zdz*vkarmn*grav/(roadrw)**1.5_wp*zalpha_w*zdL/(zus_a*zus_a*zus_a) 191 192 !! Stability function Phi_t(-z/L) (zdL is -z/L) : 193 flg = 0.5_wp + SIGN(0.5_wp, zdL) ! zdl > 0. => 1. / zdl < 0. => 0. 194 zdL2 = zdL*zdL 195 ZPHI = flg * ( 1._wp + (5._wp*zdL + 4._wp*zdL2)/(1._wp + 3._wp*zdL + 0.25_wp*zdL2) ) & ! (zdL > 0) Takaya et al. 196 & + (flg+1._wp) * ( 1._wp/SQRT(1._wp - 16._wp*(-ABS(zdL))) ) ! (zdl < 0) Eq. 8.136 197 !! FOR zdL > 0.0, old relations: 198 ! ZPHI = 1.+5._wp*zdL ! Eq. 8.136 (Large et al. 1994) 199 ! ZPHI = 1.+5.0*(zdL+zdL**2)/(1.0+3.0*zdL+zdL**2) ! SHEBA, Grachev et al. 2007 200 201 !! *** 2nd rhs term in eq. 8.156 (IFS doc Cy45r1): 202 ZL2 = - (rNu0 + 1._wp) * vkarmn * zusw / ( zdz * ZPHI ) 203 204 ! Forward time / explicit solving of eq. 8.156 (IFS doc Cy45r1): (f_n+1 == pdT(ji,jj) ; f_n == dT_wl) 205 dT_wl = MAX ( dT_wl + rdt*ZL1 + rdt*ZL2*dT_wl , 0._wp ) 206 207 ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz) 208 !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced! 209 210 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz ) ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$ 211 pdT(ji,jj) = dT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 212 183 ! Langmuir: 184 IF ( l_pustk_known ) THEN 185 zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6)) 186 ELSE 187 zla = 0.3_wp 188 END IF 189 zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp ) ! Eq.(6) 190 191 zwf = 0.5_wp + SIGN(0.5_wp, zQabs) ! zQabs > 0. => 1. / zQabs < 0. => 0. 192 193 zcst1 = vkarmn*grav*zalpha 194 195 ! 1/L when zQabs > 0 : 196 zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw) 197 198 zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 ) !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ??? 199 200 zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl 201 202 ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w ) 203 204 zcst3 = -zcst0 * vkarmn * zusw * zfLa 205 206 !! Sorry about all these constants ( constant w.r.t zdTwl), it's for 207 !! the sake of optimizations... So all these operations are not done 208 !! over and over within the iteration loop... 209 210 !! T R U L L Y I M P L I C I T => needs itteration 211 !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0.. 212 !! (without this term otherwize the implicit analytical solution is straightforward...) 213 zdTwl_n = zdTwl_b 214 DO jc = 1, 10 215 216 zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence 217 218 ! 1/L when zdTwl > 0 .AND. zQabs < 0 : 219 zL1 = SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw 220 !zL1 = vkarmn*SQRT( zdTwl_n *grav*zalpha / ( 5._wp*zHwl ) ) / zusw ! => vkarmn outside, not inside zcst1 (just for this particular line) ??? 221 222 ! Stability parameter (z/L): 223 zeta = (1._wp - zwf) * zHwl*zL1 + zwf * zHwl*zL2 224 225 ZB = zcst3 / PHI(zeta) 226 227 zdTwl_n = MAX ( zdTwl_b + ZA + ZB*zdTwl_n , 0._wp ) ! Eq.(6) 228 229 END DO 230 231 !! Update: 232 dT_wl(ji,jj) = zdTwl_n * ztcorr 233 213 234 END DO 214 235 END DO … … 216 237 END SUBROUTINE WL_ECMWF 217 238 239 240 241 FUNCTION delta_skin_layer( palpha, pQabs, pustar_a ) 242 !!--------------------------------------------------------------------- 243 !! Computes the thickness (m) of the viscous skin layer. 244 !! Based on Fairall et al., 1996 245 !! 246 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 247 !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 248 !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 249 !! doi:10.1029/95JC03190. 250 !! 251 !! L. Brodeau, october 2019 252 !!--------------------------------------------------------------------- 253 REAL(wp) :: delta_skin_layer 254 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 255 REAL(wp), INTENT(in) :: pQabs ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 256 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] 257 !!--------------------------------------------------------------------- 258 REAL(wp) :: zusw, zusw2, zlamb, zQb 259 !!--------------------------------------------------------------------- 260 261 zQb = pQabs 262 263 !zQ = MIN( -0.1_wp , pQabs ) 264 265 !ztf = 0.5_wp + SIGN(0.5_wp, zQ) ! Qabs < 0 => cooling of the layer => ztf = 0 (normal case) 266 ! ! Qabs > 0 => warming of the layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) 267 zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water 268 zusw2 = zusw*zusw 269 270 zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 271 !zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQ, 0._wp)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 272 273 delta_skin_layer = zlamb*rnu0_w/zusw 274 275 !delta_skin_layer = (1._wp - ztf) * zlamb*rnu0_w/zusw & ! see eq.(12) in Fairall et al., 1996 276 ! & + ztf * MIN(6._wp*rnu0_w/zusw , 0.007_wp) 277 END FUNCTION delta_skin_layer 278 279 280 281 FUNCTION PHI( pzeta) 282 !!--------------------------------------------------------------------- 283 !! 284 !! Takaya et al., 2010 285 !! Eq.(5) 286 !! L. Brodeau, october 2019 287 !!--------------------------------------------------------------------- 288 REAL(wp) :: PHI 289 REAL(wp), INTENT(in) :: pzeta ! stability parameter 290 !!--------------------------------------------------------------------- 291 REAL(wp) :: ztf, zzt2 292 !!--------------------------------------------------------------------- 293 ! 294 zzt2 = pzeta*pzeta 295 ! 296 ztf = 0.5_wp + SIGN(0.5_wp, pzeta) ! zeta > 0 => ztf = 1 297 ! ! zeta < 0 => ztf = 0 298 PHI = ztf * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) & ! zeta > 0 299 & + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) ) ! zeta < 0 300 ! 301 END FUNCTION PHI 302 218 303 !!====================================================================== 219 304 END MODULE sbcblk_skin_ecmwf
Note: See TracChangeset
for help on using the changeset viewer.
