Changeset 13017
 Timestamp:
 20200603T12:55:26+02:00 (16 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/branches/2020/r4.0HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk_algo_ncar.F90
r10190 r13017 11 11 !! 12 12 !! Routine turb_ncar maintained and developed in AeroBulk 13 !! (http ://aerobulk.sourceforge.net/)13 !! (https://github.com/brodeau/aerobulk/) 14 14 !! 15 !! L. Brodeau, 20 1515 !! L. Brodeau, 2020 16 16 !!===================================================================== 17 !! History : 3.6 ! 201602(L.Brodeau) successor of old turb_ncar of former sbcblk_core.F9017 !! History : 4.0 ! 202006 (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90 18 18 !! 19 19 … … 42 42 PRIVATE 43 43 44 PUBLIC :: 44 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 45 45 46 46 ! ! NCAR own values for given constants: 47 47 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature... 48 48 49 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 50 49 51 !! 50 52 CONTAINS 51 53 52 54 SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 53 & Cd, Ch, Ce, t_zu, q_zu, U _blk,&54 & Cd n, Chn, Cen)55 !! 55 & Cd, Ch, Ce, t_zu, q_zu, Ub, & 56 & CdN, ChN, CeN ) 57 !! 56 58 !! *** ROUTINE turb_ncar *** 57 59 !! … … 59 61 !! fluxes according to Large & Yeager (2004) and Large & Yeager (2008) 60 62 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 61 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 62 !! 63 !! ** Method : Monin Obukhov Similarity Theory 64 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 65 !! 66 !! ** References : Large & Yeager, 2004 / Large & Yeager, 2008 67 !! 68 !! ** Last update: Laurent Brodeau, June 2014: 69 !!  handles both cases zt=zu and zt/=zu 70 !!  optimized: less 2D arrays allocated and less operations 71 !!  better first guess of stability by checking airsea difference of virtual temperature 72 !! rather than temperature difference only... 73 !!  added function "cd_neutral_10m" that uses the improved parametrization of 74 !! Large & Yeager 2008. Dragcoefficient reduction for Cyclone conditions! 75 !!  using codewide physical constants defined into "phycst.mod" rather than redifining them 76 !! => 'vkarmn' and 'grav' 77 !! 78 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 63 !! Returns the effective bulk wind speed at zu to be used in the bulk formulas 79 64 !! 80 65 !! INPUT : 81 66 !!  82 67 !! * zt : height for temperature and spec. hum. of air [m] 83 !! * zu : height for wind speed (generally 10m) [m] 84 !! * U_zu : scalar wind speed at 10m [m/s] 85 !! * sst : SST [K] 68 !! * zu : height for wind speed (usually 10m) [m] 69 !! * sst : bulk SST [K] 86 70 !! * t_zt : potential air temperature at zt [K] 87 71 !! * ssq : specific humidity at saturation at SST [kg/kg] 88 72 !! * q_zt : specific humidity of air at zt [kg/kg] 89 !! 73 !! * U_zu : scalar wind speed at zu [m/s] 90 74 !! 91 75 !! OUTPUT : … … 96 80 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 97 81 !! * q_zu : specific humidity of air // [kg/kg] 98 !! * U_blk : bulk wind at 10m [m/s] 82 !! * Ub : bulk wind speed at zu [m/s] 83 !! 84 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 99 85 !! 100 86 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 103 89 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 90 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq ! sea surface specific humidity [kg/kg] 105 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity 91 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 92 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 107 93 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) … … 110 96 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 97 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] 112 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind at 10m [m/s] 113 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 ! 115 INTEGER :: j_itt 116 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 117 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 98 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ub ! bulk wind speed at zu [m/s] 99 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: CdN, ChN, CeN ! neutral transfer coefficients 100 ! 101 INTEGER :: j_itt 102 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 118 103 ! 119 104 REAL(wp), DIMENSION(jpi,jpj) :: Cx_n10 ! 10m neutral latent/sensible coefficient 120 REAL(wp), DIMENSION(jpi,jpj) :: sqrt _Cd_n10 ! root square of Cd_n10105 REAL(wp), DIMENSION(jpi,jpj) :: sqrtCdN10 ! root square of Cd_n10 121 106 REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu 122 107 REAL(wp), DIMENSION(jpi,jpj) :: zpsi_h_u 123 108 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 124 REAL(wp), DIMENSION(jpi,jpj) :: stab ! stability test integer 125 !! 126 ! 127 l_zt_equal_zu = .FALSE. 128 IF( ABS(zu  zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 129 130 U_blk = MAX( 0.5 , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 131 132 !! First guess of stability: 133 ztmp0 = t_zt*(1. + rctv0*q_zt)  sst*(1. + rctv0*ssq) ! airsea difference of virtual pot. temp. at zt 134 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable 135 136 !! Neutral coefficients at 10m: 109 REAL(wp), DIMENSION(jpi,jpj) :: sqrtCd 110 !! 111 112 l_zt_equal_zu = ( ABS(zu  zt) < 0.01_wp ) 113 114 Ub = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 115 116 !! Neutral drag coefficient at zu: 137 117 IF( ln_cdgw ) THEN ! wave drag case 138 cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp  tmask(:,:,1) ) 139 ztmp0 (:,:) = cdn_wave(:,:) 118 CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp  tmask(:,:,1) ) , 0.1E3_wp ) 140 119 ELSE 141 ztmp0 = cd_neutral_10m( U_blk)120 CdN = CD_N10_NCAR( Ub ) 142 121 ENDIF 143 144 sqrt_Cd_n10 = SQRT( ztmp0 ) 122 sqrtCdN10 = SQRT( CdN ) 145 123 146 124 !! Initializing transf. coeff. with their first guess neutral equivalents : 147 Cd = ztmp0 148 Ce = 1.e3*( 34.6 * sqrt_Cd_n10 ) 149 Ch = 1.e3*sqrt_Cd_n10*(18.*stab + 32.7*(1.  stab)) 150 stab = sqrt_Cd_n10 ! Temporaty array !!! stab == SQRT(Cd) 151 152 IF( ln_cdgw ) Cen = Ce ; Chn = Ch 125 Cd = CdN 126 Ce = CE_N10_NCAR( sqrtCdN10 ) 127 ztmp0 = 0.5_wp + SIGN(0.5_wp, virt_temp(t_zt, q_zt)  virt_temp(sst, ssq)) ! we guess stability based on delta of virt. pot. temp. 128 Ch = CH_N10_NCAR( sqrtCdN10 , ztmp0 ) 129 sqrtCd = sqrtCdN10 153 130 154 131 !! Initializing values at z_u with z_t values: 155 t_zu = t_zt ; q_zu = q_zt 156 157 !! * Now starting iteration loop 158 DO j_itt=1, nb_itt 132 t_zu = t_zt 133 q_zu = q_zt 134 135 !! ITERATION BLOCK 136 DO j_itt = 1, nb_itt 159 137 ! 160 138 ztmp1 = t_zu  sst ! Updating air/sea differences 161 139 ztmp2 = q_zu  ssq 162 140 163 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 164 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 165 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 166 167 ztmp0 = 1. + rctv0*q_zu ! multiply this with t and you have the virtual temperature 168 169 ! Estimate the inverse of MoninObukov length (1/L) at height zu: 170 ztmp0 = (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 171 ! ( Cd*U_blk*U_blk is U*^2 at zu ) 141 ! Updating turbulent scales : (L&Y 2004 Eq. (7)) 142 ztmp0 = sqrtCd*Ub ! u* 143 ztmp1 = Ch/sqrtCd*ztmp1 ! theta* 144 ztmp2 = Ce/sqrtCd*ztmp2 ! q* 145 146 ! Estimate the inverse of Obukov length (1/L) at height zu: 147 ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 172 148 173 149 !! Stability parameters : 174 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )175 z psi_h_u = psi_h(zeta_u )176 177 !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b9c))150 zeta_u = zu*ztmp0 151 zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 152 153 !! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b9c)) 178 154 IF( .NOT. l_zt_equal_zu ) THEN 179 !! Array 'stab' is free for the moment so using it to store 'zeta_t' 180 stab = zt*ztmp0 ; stab = SIGN( MIN(ABS(stab),10.0), stab ) ! Temporaty array stab == zeta_t !!! 181 stab = LOG(zt/zu) + zpsi_h_u  psi_h(stab) ! stab just used as temp array again! 182 t_zu = t_zt  ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b) 183 q_zu = q_zt  ztmp2/vkarmn*stab ! ztmp2 is still q* L&Y 2004 eq.(9c) 184 q_zu = max(0., q_zu) 155 ztmp0 = zt*ztmp0 ! zeta_t ! 156 ztmp0 = SIGN( MIN(ABS(ztmp0),10._wp), ztmp0 ) ! Temporaty array ztmp0 == zeta_t !!! 157 ztmp0 = LOG(zt/zu) + psi_h_ncar(zeta_u)  psi_h_ncar(ztmp0) ! ztmp0 just used as temp array again! 158 t_zu = t_zt  ztmp1/vkarmn*ztmp0 ! ztmp1 is still theta* L&Y 2004 Eq. (9b) 159 !! 160 q_zu = q_zt  ztmp2/vkarmn*ztmp0 ! ztmp2 is still q* L&Y 2004 Eq. (9c) 161 q_zu = MAX(0._wp, q_zu) 185 162 END IF 186 163 187 ztmp2 = psi_m(zeta_u) 188 IF( ln_cdgw ) THEN ! surface wave case 189 stab = vkarmn / ( vkarmn / sqrt_Cd_n10  ztmp2 ) ! (stab == SQRT(Cd)) 190 Cd = stab * stab 191 ztmp0 = (LOG(zu/10.)  zpsi_h_u) / vkarmn / sqrt_Cd_n10 192 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 193 ztmp1 = 1. + Chn * ztmp0 194 Ch = Chn * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 195 ztmp1 = 1. + Cen * ztmp0 196 Ce = Cen * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 197 164 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)... 165 ! In very rare lowwind conditions, the old way of estimating the 166 ! neutral wind speed at 10m leads to a negative value that causes the code 167 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 168 ztmp2 = psi_m_ncar(zeta_u) 169 ztmp0 = MAX( 0.25_wp , UN10_from_CD(zu, Ub, Cd, ppsi=ztmp2) ) ! U_n10 (ztmp2 == psi_m_ncar(zeta_u)) 170 171 IF( ln_cdgw ) THEN ! wave drag case 172 CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp  tmask(:,:,1) ) , 0.1E3_wp ) 198 173 ELSE 199 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 200 ! In very rare lowwind conditions, the old way of estimating the 201 ! neutral wind speed at 10m leads to a negative value that causes the code 202 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 203 ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.)  ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 204 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 205 Cdn(:,:) = ztmp0 206 sqrt_Cd_n10 = sqrt(ztmp0) 207 208 stab = 0.5 + sign(0.5,zeta_u) ! update stability 209 Cx_n10 = 1.e3*sqrt_Cd_n10*(18.*stab + 32.7*(1.  stab)) ! L&Y 2004 eq. (6c6d) (Cx_n10 == Ch_n10) 210 Chn(:,:) = Cx_n10 211 212 !! Update of transfer coefficients: 213 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.)  ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 214 Cd = ztmp0 / ( ztmp1*ztmp1 ) 215 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 216 217 ztmp0 = (LOG(zu/10.)  zpsi_h_u) / vkarmn / sqrt_Cd_n10 218 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 219 ztmp1 = 1. + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10) 220 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 221 222 Cx_n10 = 1.e3 * (34.6 * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 223 Cen(:,:) = Cx_n10 224 ztmp1 = 1. + Cx_n10*ztmp0 225 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 226 ENDIF 227 ! 228 END DO 229 ! 174 CdN = CD_N10_NCAR(ztmp0) ! Cd_n10 175 END IF 176 sqrtCdN10 = SQRT(CdN) 177 178 !! Update of transfer coefficients: 179 ztmp1 = 1._wp + sqrtCdN10/vkarmn*(LOG(zu/10._wp)  ztmp2) ! L&Y 2004 Eq. (10a) (ztmp2 == psi_m(zeta_u)) 180 Cd = MAX( CdN / ( ztmp1*ztmp1 ) , 0.1E3_wp ) 181 sqrtCd = SQRT( Cd ) 182 183 ztmp0 = ( LOG(zu/10._wp)  psi_h_ncar(zeta_u) ) / vkarmn / sqrtCdN10 184 ztmp2 = sqrtCd / sqrtCdN10 185 186 ztmp1 = 0.5_wp + sign(0.5_wp,zeta_u) ! stability flag 187 ChN = CH_N10_NCAR( sqrtCdN10 , ztmp1 ) 188 ztmp1 = 1._wp + ChN*ztmp0 189 Ch = MAX( ChN*ztmp2 / ztmp1 , 0.1E3_wp ) ! L&Y 2004 Eq. (10b) 190 191 CeN = CE_N10_NCAR( sqrtCdN10 ) 192 ztmp1 = 1._wp + CeN*ztmp0 193 Ce = MAX( CeN*ztmp2 / ztmp1 , 0.1E3_wp ) ! L&Y 2004 Eq. (10c) 194 195 END DO !DO j_itt = 1, nb_itt 196 230 197 END SUBROUTINE turb_ncar 231 198 232 199 233 FUNCTION cd_neutral_10m( pw10 )234 !! 200 FUNCTION CD_N10_NCAR( pw10 ) 201 !! 235 202 !! Estimate of the neutral drag coefficient at 10m as a function 236 203 !! of neutral wind speed at 10m 237 204 !! 238 !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)239 !! 240 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)205 !! Origin: Large & Yeager 2008, Eq. (11) 206 !! 207 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 241 208 !! 242 209 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s) 243 REAL(wp), DIMENSION(jpi,jpj) :: cd_neutral_10m210 REAL(wp), DIMENSION(jpi,jpj) :: CD_N10_NCAR 244 211 ! 245 212 INTEGER :: ji, jj ! dummy loop indices … … 255 222 ! 256 223 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 257 zgt33 = 0.5 + SIGN( 0.5, (zw  33.) ) ! If pw10 < 33. => 0, else => 1258 ! 259 cd_neutral_10m(ji,jj) = 1.e3* ( &260 & (1.  zgt33)*( 2.7/zw + 0.142 + zw/13.09  3.14807E10*zw6) & ! wind < 33 m/s261 & + zgt33 * 2.34 )! wind >= 33 m/s262 ! 263 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E6)224 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw  33._wp) ) ! If pw10 < 33. => 0, else => 1 225 ! 226 CD_N10_NCAR(ji,jj) = 1.e3_wp * ( & 227 & (1._wp  zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp  3.14807E10_wp*zw6) & ! wind < 33 m/s 228 & + zgt33 * 2.34_wp ) ! wind >= 33 m/s 229 ! 230 CD_N10_NCAR(ji,jj) = MAX( CD_N10_NCAR(ji,jj), 0.1E3_wp ) 264 231 ! 265 232 END DO 266 233 END DO 267 234 ! 268 END FUNCTION cd_neutral_10m 269 270 271 FUNCTION psi_m( pzeta ) 235 END FUNCTION CD_N10_NCAR 236 237 238 239 FUNCTION CH_N10_NCAR( psqrtcdn10 , pstab ) 240 !! 241 !! Estimate of the neutral heat transfer coefficient at 10m !! 242 !! Origin: Large & Yeager 2008, Eq. (9) and (12) 243 244 !! 245 REAL(wp), DIMENSION(jpi,jpj) :: ch_n10_ncar 246 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 247 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab ! stable ABL => 1 / unstable ABL => 0 248 !! 249 ! 250 ch_n10_ncar = MAX( 1.e3_wp * psqrtcdn10*( 18._wp*pstab + 32.7_wp*(1._wp  pstab) ) , 0.1E3_wp ) ! Eq. (9) & (12) Large & Yeager, 2008 251 ! 252 END FUNCTION CH_N10_NCAR 253 254 FUNCTION CE_N10_NCAR( psqrtcdn10 ) 255 !! 256 !! Estimate of the neutral heat transfer coefficient at 10m !! 257 !! Origin: Large & Yeager 2008, Eq. (9) and (13) 258 !! 259 REAL(wp), DIMENSION(jpi,jpj) :: ce_n10_ncar 260 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 261 !! 262 ce_n10_ncar = MAX( 1.e3_wp * ( 34.6_wp * psqrtcdn10 ) , 0.1E3_wp ) 263 ! 264 END FUNCTION CE_N10_NCAR 265 266 267 FUNCTION psi_m_ncar( pzeta ) 272 268 !! 273 269 !! Universal profile stability function for momentum 274 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e)275 !! 276 !! pzet 0 : stability paramenter, z/L where z is altitude measurement270 !! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 271 !! 272 !! pzeta : stability paramenter, z/L where z is altitude measurement 277 273 !! and L is MO length 278 274 !! 279 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 280 !! 281 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 282 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 283 ! 284 INTEGER :: ji, jj ! dummy loop indices 285 REAL(wp) :: zx2, zx, zstab ! local scalars 286 !! 287 ! 275 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 276 !! 277 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar 278 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 279 ! 280 INTEGER :: ji, jj ! dummy loop indices 281 REAL(wp) :: zzeta, zx2, zx, zpsi_unst, zpsi_stab, zstab ! local scalars 282 !! 288 283 DO jj = 1, jpj 289 284 DO ji = 1, jpi 290 zx2 = SQRT( ABS( 1.  16.*pzeta(ji,jj) ) ) 291 zx2 = MAX ( zx2 , 1. ) 292 zx = SQRT( zx2 ) 293 zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 294 ! 295 psi_m(ji,jj) = zstab * (5.*pzeta(ji,jj)) & ! Stable 296 & + (1.  zstab) * (2.*LOG((1. + zx)*0.5) & ! Unstable 297 & + LOG((1. + zx2)*0.5)  2.*ATAN(zx) + rpi*0.5) ! " 285 286 zzeta = pzeta(ji,jj) 287 ! 288 zx2 = SQRT( ABS(1._wp  16._wp*zzeta) ) ! (1  16z)^0.5 289 zx2 = MAX( zx2 , 1._wp ) 290 zx = SQRT(zx2) ! (1  16z)^0.25 291 zpsi_unst = 2._wp*LOG( (1._wp + zx )*0.5_wp ) & 292 & + LOG( (1._wp + zx2)*0.5_wp ) & 293 &  2._wp*ATAN(zx) + rpi*0.5_wp 294 ! 295 zpsi_stab = 5._wp*zzeta 296 ! 297 zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1 298 ! 299 psi_m_ncar(ji,jj) = zstab * zpsi_stab & ! (zzeta > 0) Stable 300 & + (1._wp  zstab) * zpsi_unst ! (zzeta < 0) Unstable 298 301 ! 299 302 END DO 300 303 END DO 301 ! 302 END FUNCTION psi_m 303 304 305 FUNCTION psi_h( pzeta ) 304 END FUNCTION psi_m_ncar 305 306 307 FUNCTION psi_h_ncar( pzeta ) 306 308 !! 307 309 !! Universal profile stability function for temperature and humidity 308 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e)309 !! 310 !! pzet 0 : stability paramenter, z/L where z is altitude measurement310 !! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 311 !! 312 !! pzeta : stability paramenter, z/L where z is altitude measurement 311 313 !! and L is MO length 312 314 !! 313 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 314 !! 315 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 316 !! 317 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar 315 318 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 316 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 317 ! 318 INTEGER :: ji, jj ! dummy loop indices 319 REAL(wp) :: zx2, zstab ! local scalars 319 ! 320 INTEGER :: ji, jj ! dummy loop indices 321 REAL(wp) :: zzeta, zx2, zpsi_unst, zpsi_stab, zstab ! local scalars 320 322 !! 321 323 ! 322 324 DO jj = 1, jpj 323 325 DO ji = 1, jpi 324 zx2 = SQRT( ABS( 1.  16.*pzeta(ji,jj) ) ) 325 zx2 = MAX ( zx2 , 1. ) 326 zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 327 ! 328 psi_h(ji,jj) = zstab * (5.*pzeta(ji,jj)) & ! Stable 329 & + (1.  zstab) * (2.*LOG( (1. + zx2)*0.5 )) ! Unstable 326 ! 327 zzeta = pzeta(ji,jj) 328 ! 329 zx2 = SQRT( ABS(1._wp  16._wp*zzeta) ) ! (1 16z)^0.5 330 zx2 = MAX( zx2 , 1._wp ) 331 zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) 332 ! 333 zpsi_stab = 5._wp*zzeta 334 ! 335 zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1 336 ! 337 psi_h_ncar(ji,jj) = zstab * zpsi_stab & ! (zzeta > 0) Stable 338 & + (1._wp  zstab) * zpsi_unst ! (zzeta < 0) Unstable 330 339 ! 331 340 END DO 332 341 END DO 333 ! 334 END FUNCTION psi_h 342 END FUNCTION psi_h_ncar 343 344 345 346 347 FUNCTION UN10_from_CD( pzu, pUb, pCd, ppsi ) 348 !! 349 !! Provides the neutralstability wind speed at 10 m 350 !! 351 REAL(wp), DIMENSION(jpi,jpj) :: UN10_from_CD !: [m/s] 352 REAL(wp), INTENT(in) :: pzu !: measurement heigh of bulk wind speed 353 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb !: bulk wind speed at height pzu m [m/s] 354 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd !: drag coefficient 355 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum [] 356 !! 357 !! Reminder: UN10 = u*/vkarmn * log(10/z0) 358 !! and: u* = sqrt(Cd) * Ub 359 !! u*/vkarmn * log( 10 / z0 ) 360 UN10_from_CD(:,:) = SQRT(pCd(:,:))*pUb/vkarmn * LOG( 10._wp / z0_from_Cd( pzu, pCd(:,:), ppsi=ppsi(:,:) ) ) 361 !! 362 END FUNCTION UN10_from_CD 363 364 365 FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) 366 !! 367 !! 368 !! Evaluates the 1./(Obukhov length) from air temperature, 369 !! air specific humidity, and frictional scales u*, t* and q* 370 !! 371 !! Author: L. Brodeau, June 2019 / AeroBulk 372 !! (https://github.com/brodeau/aerobulk/) 373 !! 374 REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Obukhov length) [m^1] 375 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha !: reference potential temperature of air [K] 376 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: reference specific humidity of air [kg/kg] 377 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus !: u*: friction velocity [m/s] 378 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts, pqs !: \theta* and q* friction aka turb. scales for temp. and spec. hum. 379 ! 380 INTEGER :: ji, jj ! dummy loop indices 381 REAL(wp) :: zqa ! local scalar 382 !! 383 ! 384 DO jj = 1, jpj 385 DO ji = 1, jpi 386 ! 387 zqa = (1._wp + rctv0*pqa(ji,jj)) 388 ! 389 ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 390 ! a/ u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 391 ! or 392 ! b/ u* [ theta* + 0.61 theta q* ] 393 ! 394 One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 395 & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E9_wp ) 396 ! 397 END DO 398 END DO 399 ! 400 One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...) 401 ! 402 END FUNCTION One_on_L 403 404 405 FUNCTION z0_from_Cd( pzu, pCd, ppsi ) 406 REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd !: roughness length [m] 407 REAL(wp) , INTENT(in) :: pzu !: reference height zu [m] 408 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd !: (neutral or nonneutral) drag coefficient [] 409 REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum [] 410 !! 411 !! If pCd is the NEUTRALSTABILITY drag coefficient then ppsi must be 0 or not given 412 !! If pCd is the drag coefficient (in stable or unstable conditions) then pssi must be provided 413 !! 414 IF ( PRESENT(ppsi) ) THEN 415 !! Cd provided is the actual Cd (not the neutralstability CdN) : 416 z0_from_Cd = pzu * EXP(  ( vkarmn/SQRT(pCd(:,:)) + ppsi(:,:) ) ) !LB: ok, doublechecked! 417 ELSE 418 !! Cd provided is the neutralstability Cd, aka CdN : 419 z0_from_Cd = pzu * EXP(  vkarmn/SQRT(pCd(:,:)) ) !LB: ok, doublechecked! 420 END IF 421 END FUNCTION z0_from_Cd 422 423 FUNCTION virt_temp( pta, pqa ) 424 REAL(wp), DIMENSION(jpi,jpj) :: virt_temp !: virtual temperature [K] 425 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K] 426 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air [kg/kg] 427 virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) 428 END FUNCTION virt_temp 335 429 336 430 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.