Changeset 11209 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
- Timestamp:
- 2019-07-03T09:09:17+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
r11111 r11209 38 38 USE lib_fortran ! to use key_nosignedzero 39 39 40 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 40 41 41 42 IMPLICIT NONE … … 93 94 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 95 !! * q_zu : specific humidity of air // [kg/kg] 95 !! * U_blk : bulk wind at 10m[m/s]96 !! * U_blk : bulk wind speed at 10m [m/s] 96 97 !!---------------------------------------------------------------------------------- 97 98 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 125 126 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 126 127 127 U_blk = MAX( 0.5 , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s128 U_blk = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 128 129 129 130 !! First guess of stability: 130 ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt131 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable131 ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 132 stab = 0.5_wp + sign(0.5_wp,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable 132 133 133 134 !! Neutral coefficients at 10m: … … 159 160 160 161 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 161 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 162 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 163 164 ztmp0 = 1. + rctv0*q_zu ! multiply this with t and you have the virtual temperature 162 ztmp0 = stab*U_blk ! u* (stab == SQRT(Cd)) 163 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 164 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 165 165 166 166 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 167 ztmp0 = (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 168 ! ( Cd*U_blk*U_blk is U*^2 at zu ) 169 167 ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 168 170 169 !! Stability parameters : 171 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 170 zeta_u = zu*ztmp0 171 zeta_u = sign( min(abs(zeta_u),10.0_wp), zeta_u ) 172 172 zpsi_h_u = psi_h( zeta_u ) 173 173 … … 175 175 IF( .NOT. l_zt_equal_zu ) THEN 176 176 !! Array 'stab' is free for the moment so using it to store 'zeta_t' 177 stab = zt*ztmp0 ; stab = SIGN( MIN(ABS(stab),10.0), stab ) ! Temporaty array stab == zeta_t !!! 177 stab = zt*ztmp0 178 stab = SIGN( MIN(ABS(stab),10.0_wp), stab ) ! Temporaty array stab == zeta_t !!! 178 179 stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again! 179 180 t_zu = t_zt - ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b) 180 181 q_zu = q_zt - ztmp2/vkarmn*stab ! ztmp2 is still q* L&Y 2004 eq.(9c) 181 q_zu = max(0. , q_zu)182 q_zu = max(0._wp, q_zu) 182 183 END IF 183 184 … … 194 195 195 196 ELSE 196 197 198 199 200 ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))201 202 203 204 205 stab = 0.5 + sign(0.5,zeta_u)! update stability206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 197 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 198 ! In very rare low-wind conditions, the old way of estimating the 199 ! neutral wind speed at 10m leads to a negative value that causes the code 200 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 201 ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 202 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 203 Cdn(:,:) = ztmp0 204 sqrt_Cd_n10 = sqrt(ztmp0) 205 206 stab = 0.5_wp + sign(0.5_wp,zeta_u) ! update stability 207 Cx_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10) 208 Chn(:,:) = Cx_n10 209 210 !! Update of transfer coefficients: 211 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 212 Cd = ztmp0 / ( ztmp1*ztmp1 ) 213 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 214 215 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 216 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 217 ztmp1 = 1. + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10) 218 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 219 220 Cx_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 221 Cen(:,:) = Cx_n10 222 ztmp1 = 1. + Cx_n10*ztmp0 223 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 224 ENDIF 224 225 ! 225 226 END DO … … 252 253 ! 253 254 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 254 zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) ) ! If pw10 < 33. => 0, else => 1255 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1 255 256 ! 256 257 cd_neutral_10m(ji,jj) = 1.e-3 * ( & … … 258 259 & + zgt33 * 2.34 ) ! wind >= 33 m/s 259 260 ! 260 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6 )261 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 261 262 ! 262 263 END DO … … 285 286 DO jj = 1, jpj 286 287 DO ji = 1, jpi 287 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )288 zx2 = MAX ( zx2 , 1.)288 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 289 zx2 = MAX( zx2 , 1._wp ) 289 290 zx = SQRT( zx2 ) 290 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )291 ! 292 psi_m(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable293 & + (1. - zstab) * (2.*LOG((1. + zx)*0.5) & ! Unstable294 & + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5) ! "291 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 292 ! 293 psi_m(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 294 & + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp) & ! Unstable 295 & + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp) ! " 295 296 ! 296 297 END DO … … 319 320 DO jj = 1, jpj 320 321 DO ji = 1, jpi 321 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )322 zx2 = MAX ( zx2 , 1.)323 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )324 ! 325 psi_h(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable326 & + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5)) ! Unstable322 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 323 zx2 = MAX( zx2 , 1._wp ) 324 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 325 ! 326 psi_h(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 327 & + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp )) ! Unstable 327 328 ! 328 329 END DO
Note: See TracChangeset
for help on using the changeset viewer.
