Changeset 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk_algo_ncar.F90
- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/SBC/sbcblk_algo_ncar.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk_algo_ncar.F90
r10190 r12928 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 15 !! L. Brodeau, 2015 … … 38 38 USE lib_fortran ! to use key_nosignedzero 39 39 40 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 40 41 41 42 IMPLICIT NONE 42 43 PRIVATE 43 44 44 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 45 46 ! ! NCAR own values for given constants: 47 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature... 48 45 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 46 47 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 48 !! * Substitutions 49 # include "do_loop_substitute.h90" 50 49 51 !!---------------------------------------------------------------------- 50 52 CONTAINS … … 61 63 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 62 64 !! 63 !! ** Method : Monin Obukhov Similarity Theory64 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)65 !!66 !! ** References : Large & Yeager, 2004 / Large & Yeager, 200867 !!68 !! ** Last update: Laurent Brodeau, June 2014:69 !! - handles both cases zt=zu and zt/=zu70 !! - optimized: less 2D arrays allocated and less operations71 !! - better first guess of stability by checking air-sea difference of virtual temperature72 !! rather than temperature difference only...73 !! - added function "cd_neutral_10m" that uses the improved parametrization of74 !! Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!75 !! - using code-wide physical constants defined into "phycst.mod" rather than redifining them76 !! => 'vkarmn' and 'grav'77 !!78 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)79 65 !! 80 66 !! INPUT : 81 67 !! ------- 82 68 !! * 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] 69 !! * zu : height for wind speed (usually 10m) [m] 70 !! * sst : bulk SST [K] 86 71 !! * t_zt : potential air temperature at zt [K] 87 72 !! * ssq : specific humidity at saturation at SST [kg/kg] 88 73 !! * q_zt : specific humidity of air at zt [kg/kg] 74 !! * U_zu : scalar wind speed at zu [m/s] 89 75 !! 90 76 !! … … 96 82 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 97 83 !! * q_zu : specific humidity of air // [kg/kg] 98 !! * U_blk : bulk wind at 10m [m/s] 84 !! * U_blk : bulk wind speed at zu [m/s] 85 !! 86 !! 87 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 99 88 !!---------------------------------------------------------------------------------- 100 89 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 103 92 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 93 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 94 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 95 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 107 96 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) … … 110 99 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 100 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]101 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 113 102 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 103 ! 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 104 INTEGER :: j_itt 105 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 118 106 ! 119 107 REAL(wp), DIMENSION(jpi,jpj) :: Cx_n10 ! 10m neutral latent/sensible coefficient … … 124 112 REAL(wp), DIMENSION(jpi,jpj) :: stab ! stability test integer 125 113 !!---------------------------------------------------------------------------------- 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 114 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 115 116 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 131 117 132 118 !! First guess of stability: 133 ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt134 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable119 ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 120 stab = 0.5_wp + sign(0.5_wp,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable 135 121 136 122 !! Neutral coefficients at 10m: … … 139 125 ztmp0 (:,:) = cdn_wave(:,:) 140 126 ELSE 141 127 ztmp0 = cd_neutral_10m( U_blk ) 142 128 ENDIF 143 129 … … 146 132 !! Initializing transf. coeff. with their first guess neutral equivalents : 147 133 Cd = ztmp0 148 Ce = 1.e-3 *( 34.6* sqrt_Cd_n10 )149 Ch = 1.e-3 *sqrt_Cd_n10*(18.*stab + 32.7*(1.- stab))134 Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 ) 135 Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) 150 136 stab = sqrt_Cd_n10 ! Temporaty array !!! stab == SQRT(Cd) 151 137 152 IF( ln_cdgw ) Cen = Ce ; Chn = Ch 153 154 !! 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 138 IF( ln_cdgw ) THEN 139 Cen = Ce 140 Chn = Ch 141 ENDIF 142 143 !! First guess of temperature and humidity at height zu: 144 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 145 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 146 147 !! ITERATION BLOCK 148 DO j_itt = 1, nb_itt 159 149 ! 160 150 ztmp1 = t_zu - sst ! Updating air/sea differences … … 162 152 163 153 ! 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 154 ztmp0 = stab*U_blk ! u* (stab == SQRT(Cd)) 155 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 156 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 168 157 169 158 ! Estimate the inverse of Monin-Obukov 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 ) 172 159 ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 160 173 161 !! Stability parameters : 174 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 162 zeta_u = zu*ztmp0 163 zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 175 164 zpsi_h_u = psi_h( zeta_u ) 176 165 … … 178 167 IF( .NOT. l_zt_equal_zu ) THEN 179 168 !! 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 !!! 169 stab = zt*ztmp0 170 stab = SIGN( MIN(ABS(stab),10._wp), stab ) ! Temporaty array stab == zeta_t !!! 181 171 stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again! 182 172 t_zu = t_zt - ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b) 183 173 q_zu = q_zt - ztmp2/vkarmn*stab ! ztmp2 is still q* L&Y 2004 eq.(9c) 184 q_zu = max(0., q_zu) 185 END IF 186 174 q_zu = max(0._wp, q_zu) 175 ENDIF 176 177 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 178 ! In very rare low-wind conditions, the old way of estimating the 179 ! neutral wind speed at 10m leads to a negative value that causes the code 180 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 187 181 ztmp2 = psi_m(zeta_u) 188 182 IF( ln_cdgw ) THEN ! surface wave case 189 183 stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 ) ! (stab == SQRT(Cd)) 190 184 Cd = stab * stab 191 ztmp0 = (LOG(zu/10. ) - zpsi_h_u) / vkarmn / sqrt_Cd_n10185 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 192 186 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 193 ztmp1 = 1. + Chn * ztmp0187 ztmp1 = 1._wp + Chn * ztmp0 194 188 Ch = Chn * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 195 ztmp1 = 1. + Cen * ztmp0189 ztmp1 = 1._wp + Cen * ztmp0 196 190 Ce = Cen * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 197 191 198 192 ELSE 199 200 201 202 203 ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))204 205 206 207 208 stab = 0.5 + sign(0.5,zeta_u)! update stability209 Cx_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.- stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10)210 211 212 213 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))214 215 216 217 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10218 219 ztmp1 = 1.+ Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10)220 221 222 Cx_n10 = 1.e-3 * (34.6* sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10223 224 ztmp1 = 1.+ Cx_n10*ztmp0225 226 227 ! 228 END DO 229 ! 193 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 194 ! In very rare low-wind conditions, the old way of estimating the 195 ! neutral wind speed at 10m leads to a negative value that causes the code 196 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 197 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)) 198 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 199 Cdn(:,:) = ztmp0 200 sqrt_Cd_n10 = sqrt(ztmp0) 201 202 stab = 0.5_wp + sign(0.5_wp,zeta_u) ! update stability 203 Cx_n10 = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10) 204 Chn(:,:) = Cx_n10 205 206 !! Update of transfer coefficients: 207 ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 208 Cd = ztmp0 / ( ztmp1*ztmp1 ) 209 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 210 211 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 212 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 213 ztmp1 = 1._wp + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10) 214 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 215 216 Cx_n10 = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 217 Cen(:,:) = Cx_n10 218 ztmp1 = 1._wp + Cx_n10*ztmp0 219 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 220 ENDIF 221 222 END DO !DO j_itt = 1, nb_itt 223 230 224 END SUBROUTINE turb_ncar 231 225 … … 238 232 !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 239 233 !! 240 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)234 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 241 235 !!---------------------------------------------------------------------------------- 242 236 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s) … … 247 241 !!---------------------------------------------------------------------------------- 248 242 ! 249 DO jj = 1, jpj 250 DO ji = 1, jpi 251 ! 252 zw = pw10(ji,jj) 253 zw6 = zw*zw*zw 254 zw6 = zw6*zw6 255 ! 256 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 257 zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) ) ! If pw10 < 33. => 0, else => 1 258 ! 259 cd_neutral_10m(ji,jj) = 1.e-3 * ( & 260 & (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind < 33 m/s 261 & + zgt33 * 2.34 ) ! wind >= 33 m/s 262 ! 263 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6) 264 ! 265 END DO 266 END DO 243 DO_2D_11_11 244 ! 245 zw = pw10(ji,jj) 246 zw6 = zw*zw*zw 247 zw6 = zw6*zw6 248 ! 249 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 250 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1 251 ! 252 cd_neutral_10m(ji,jj) = 1.e-3_wp * ( & 253 & (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind < 33 m/s 254 & + zgt33 * 2.34_wp ) ! wind >= 33 m/s 255 ! 256 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 257 ! 258 END_2D 267 259 ! 268 260 END FUNCTION cd_neutral_10m … … 273 265 !! Universal profile stability function for momentum 274 266 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 275 !! 276 !! pzet 0 : stability paramenter, z/L where z is altitude measurement267 !! 268 !! pzeta : stability paramenter, z/L where z is altitude measurement 277 269 !! and L is M-O length 278 270 !! 279 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)280 !!---------------------------------------------------------------------------------- 281 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pzeta282 REAL(wp), DIMENSION(jpi,jpj) :: psi_m283 ! 284 INTEGER :: ji, jj 271 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 272 !!---------------------------------------------------------------------------------- 273 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 274 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 275 ! 276 INTEGER :: ji, jj ! dummy loop indices 285 277 REAL(wp) :: zx2, zx, zstab ! local scalars 286 278 !!---------------------------------------------------------------------------------- 287 ! 288 DO jj = 1, jpj 289 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) ! " 298 ! 299 END DO 300 END DO 301 ! 279 DO_2D_11_11 280 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 281 zx2 = MAX( zx2 , 1._wp ) 282 zx = SQRT( zx2 ) 283 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 284 ! 285 psi_m(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 286 & + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp) & ! Unstable 287 & + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp) ! " 288 ! 289 END_2D 302 290 END FUNCTION psi_m 303 291 … … 308 296 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 309 297 !! 310 !! pzet 0 : stability paramenter, z/L where z is altitude measurement298 !! pzeta : stability paramenter, z/L where z is altitude measurement 311 299 !! and L is M-O length 312 300 !! 313 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 314 !!---------------------------------------------------------------------------------- 301 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 302 !!---------------------------------------------------------------------------------- 303 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 315 304 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 316 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 317 ! 318 INTEGER :: ji, jj ! dummy loop indices 305 ! 306 INTEGER :: ji, jj ! dummy loop indices 319 307 REAL(wp) :: zx2, zstab ! local scalars 320 308 !!---------------------------------------------------------------------------------- 321 309 ! 322 DO jj = 1, jpj 323 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 330 ! 331 END DO 332 END DO 333 ! 310 DO_2D_11_11 311 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 312 zx2 = MAX( zx2 , 1._wp ) 313 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 314 ! 315 psi_h(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 316 & + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp )) ! Unstable 317 ! 318 END_2D 334 319 END FUNCTION psi_h 335 320
Note: See TracChangeset
for help on using the changeset viewer.
