Context Navigation

sbcblk_algo_ncar.F90 @ 14205

Last change on this file since 14205 was 14072, checked in by laurent, 4 years ago
Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369
Property svn:keywords set to `Id`
File size: 17.6 KB

Rev	Line
[6723]	1	MODULE sbcblk_algo_ncar
	2	!!======================================================================
	3	!! * MODULE sbcblk_algo_ncar *
	4	!! Computes:
	5	!! * bulk transfer coefficients C_D, C_E and C_H
	6	!! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
[14072]	7	!! * the effective bulk wind speed at 10m Ubzu
[6723]	8	!! => all these are used in bulk formulas in sbcblk.F90
	9	!!
	10	!! Using the bulk formulation/param. of Large & Yeager 2008
	11	!!
	12	!! Routine turb_ncar maintained and developed in AeroBulk
[12377]	13	!! (https://github.com/brodeau/aerobulk/)
[6723]	14	!!
	15	!! L. Brodeau, 2015
	16	!!=====================================================================
[6727]	17	!! History : 3.6 ! 2016-02 (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90
[14072]	18	!! 4.2 ! 2020-12 (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
[6723]	19	!!----------------------------------------------------------------------
	20
	21	!!----------------------------------------------------------------------
	22	!! turb_ncar : computes the bulk turbulent transfer coefficients
	23	!! adjusts t_air and q_air from zt to zu m
	24	!! returns the effective bulk wind speed at 10m
	25	!!----------------------------------------------------------------------
	26	USE dom_oce ! ocean space and time domain
[14072]	27	USE sbc_oce, ONLY: ln_cdgw
	28	USE sbcwave, ONLY: cdn_wave ! wave module
[6723]	29	USE phycst ! physical constants
[14072]	30	USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
[6723]	31
	32	IMPLICIT NONE
	33	PRIVATE
	34
[12377]	35	PUBLIC :: TURB_NCAR ! called by sbcblk.F90
[6723]	36
[12377]	37	!! * Substitutions
	38	# include "do_loop_substitute.h90"
	39
[6723]	40	!!----------------------------------------------------------------------
	41	CONTAINS
	42
[14072]	43	SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
	44	& Cd, Ch, Ce, t_zu, q_zu, Ubzu, &
	45	& nb_iter, CdN, ChN, CeN )
[6723]	46	!!----------------------------------------------------------------------------------
	47	!! * ROUTINE turb_ncar *
	48	!!
	49	!! ** Purpose : Computes turbulent transfert coefficients of surface
	50	!! fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
	51	!! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
[14072]	52	!! Returns the effective bulk wind speed at zu to be used in the bulk formulas
[6723]	53	!!
	54	!! INPUT :
	55	!! -------
	56	!! * zt : height for temperature and spec. hum. of air [m]
[12377]	57	!! * zu : height for wind speed (usually 10m) [m]
	58	!! * sst : bulk SST [K]
[6723]	59	!! * t_zt : potential air temperature at zt [K]
	60	!! * ssq : specific humidity at saturation at SST [kg/kg]
	61	!! * q_zt : specific humidity of air at zt [kg/kg]
[12377]	62	!! * U_zu : scalar wind speed at zu [m/s]
[6723]	63	!!
	64	!!
	65	!! OUTPUT :
	66	!! --------
	67	!! * Cd : drag coefficient
	68	!! * Ch : sensible heat coefficient
	69	!! * Ce : evaporation coefficient
	70	!! * t_zu : pot. air temperature adjusted at wind height zu [K]
	71	!! * q_zu : specific humidity of air // [kg/kg]
[14072]	72	!! * Ubzu : bulk wind speed at zu [m/s]
[12377]	73	!!
[14072]	74	!! OPTIONAL OUTPUT:
	75	!! ----------------
	76	!! * CdN : neutral-stability drag coefficient
	77	!! * ChN : neutral-stability sensible heat coefficient
	78	!! * CeN : neutral-stability evaporation coefficient
[12377]	79	!!
	80	!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	81	!!----------------------------------------------------------------------------------
	82	REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m]
	83	REAL(wp), INTENT(in ) :: zu ! height for U_zu [m]
	84	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst ! sea surface temperature [Kelvin]
	85	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin]
	86	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq ! sea surface specific humidity [kg/kg]
[12377]	87	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg]
[6723]	88	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s]
	89	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau)
	90	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens)
	91	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat)
	92	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K]
	93	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg]
[14072]	94	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s]
[6723]	95	!
[14072]	96	INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations
	97	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN
	98	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN
	99	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN
	100	!
	101	INTEGER :: nbit, jit ! iterations...
[12377]	102	LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U
[6723]	103	!
[14072]	104	REAL(wp), DIMENSION(jpi,jpj) :: zCdN, zCeN, zChN ! 10m neutral latent/sensible coefficient
	105	REAL(wp), DIMENSION(jpi,jpj) :: zsqrt_Cd, zsqrt_CdN ! root square of Cd and Cd_neutral
[9125]	106	REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
	107	REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
[6723]	108	!!----------------------------------------------------------------------------------
[14072]	109	nbit = nb_iter0
	110	IF( PRESENT(nb_iter) ) nbit = nb_iter
	111
[12615]	112	l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
[6723]	113
[14072]	114	Ubzu = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
[6723]	115
	116	!! First guess of stability:
[12377]	117	ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt
[14072]	118	ztmp1 = 0.5_wp + SIGN(0.5_wp,ztmp0) ! ztmp1 = 1 if dTv > 0 => STABLE, 0 if unstable
[6723]	119
	120	!! Neutral coefficients at 10m:
	121	IF( ln_cdgw ) THEN ! wave drag case
[7753]	122	cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
[14072]	123	zCdN (:,:) = cdn_wave(:,:)
[6723]	124	ELSE
[14072]	125	zCdN = cd_n10_ncar( Ubzu )
[6723]	126	ENDIF
	127
[14072]	128	zsqrt_CdN = SQRT( zCdN )
[6723]	129
	130	!! Initializing transf. coeff. with their first guess neutral equivalents :
[14072]	131	Cd = zCdN
	132	Ce = ce_n10_ncar( zsqrt_CdN )
	133	Ch = ch_n10_ncar( zsqrt_CdN , ztmp1 ) ! ztmp1 is stability (1/0)
	134	zsqrt_Cd = zsqrt_CdN
	135
[12377]	136	IF( ln_cdgw ) THEN
[14072]	137	zCeN = Ce
	138	zChN = Ch
[12377]	139	ENDIF
[6723]	140
[14072]	141	!! Initializing values at z_u with z_t values:
[12615]	142	t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions...
	143	q_zu = MAX( q_zt , 1.e-6_wp ) ! "
[6723]	144
[14072]	145
[12377]	146	!! ITERATION BLOCK
[14072]	147	DO jit = 1, nbit
[6723]	148	!
	149	ztmp1 = t_zu - sst ! Updating air/sea differences
	150	ztmp2 = q_zu - ssq
	151
[14072]	152	! Updating turbulent scales : (L&Y 2004 Eq. (7))
	153	ztmp0 = zsqrt_CdUbzu ! u
	154	ztmp1 = Ch/zsqrt_Cdztmp1 ! theta
	155	ztmp2 = Ce/zsqrt_Cdztmp2 ! q
[6723]	156
[14072]	157	! Estimate the inverse of Obukov length (1/L) at height zu:
[12377]	158	ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
[14072]	159
[6723]	160	!! Stability parameters :
[12377]	161	zeta_u = zu*ztmp0
[14072]	162	zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u )
[6723]	163
[14072]	164	!! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b-9c))
[6723]	165	IF( .NOT. l_zt_equal_zu ) THEN
[14072]	166	ztmp0 = zt*ztmp0 ! zeta_t !
	167	ztmp0 = SIGN( MIN(ABS(ztmp0),10._wp), ztmp0 ) ! Temporaty array ztmp0 == zeta_t !!!
	168	ztmp0 = LOG(zt/zu) + psi_h_ncar(zeta_u) - psi_h_ncar(ztmp0) ! ztmp0 just used as temp array again!
	169	t_zu = t_zt - ztmp1/vkarmnztmp0 ! ztmp1 is still theta L&Y 2004 Eq. (9b)
	170	!!
	171	q_zu = q_zt - ztmp2/vkarmnztmp0 ! ztmp2 is still q L&Y 2004 Eq. (9c)
	172	q_zu = MAX(0._wp, q_zu)
	173	END IF
[6723]	174
[14072]	175	! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)...
[12377]	176	! In very rare low-wind conditions, the old way of estimating the
	177	! neutral wind speed at 10m leads to a negative value that causes the code
	178	! to crash. To prevent this a threshold of 0.25m/s is imposed.
[14072]	179	ztmp2 = psi_m_ncar(zeta_u)
[6723]	180	IF( ln_cdgw ) THEN ! surface wave case
[14072]	181	zsqrt_Cd = vkarmn / ( vkarmn / zsqrt_CdN - ztmp2 )
	182	Cd = zsqrt_Cd * zsqrt_Cd
	183	ztmp0 = (LOG(zu/10._wp) - psi_h_ncar(zeta_u)) / vkarmn / zsqrt_CdN
	184	ztmp2 = zsqrt_Cd / zsqrt_CdN
	185	ztmp1 = 1._wp + zChN * ztmp0
	186	Ch = zChN * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b)
	187	ztmp1 = 1._wp + zCeN * ztmp0
	188	Ce = zCeN * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c)
[10190]	189
[6723]	190	ELSE
[14072]	191	ztmp0 = MAX( 0.25_wp , UN10_from_CD(zu, Ubzu, Cd, ppsi=ztmp2) ) ! U_n10 (ztmp2 == psi_m_ncar(zeta_u))
[6723]	192
[14072]	193	zCdN = cd_n10_ncar(ztmp0)
	194	zsqrt_CdN = sqrt(zCdN)
[6723]	195
[12377]	196	!! Update of transfer coefficients:
[6723]	197
[14072]	198	!! C_D
	199	ztmp1 = 1._wp + zsqrt_CdN/vkarmn*(LOG(zu/10._wp) - ztmp2) ! L&Y 2004 Eq. (10a) (ztmp2 == psi_m(zeta_u))
	200	Cd = MAX( zCdN / ( ztmp1*ztmp1 ), Cx_min )
[6723]	201
[14072]	202	!! C_H and C_E
	203	zsqrt_Cd = SQRT( Cd )
	204	ztmp0 = ( LOG(zu/10._wp) - psi_h_ncar(zeta_u) ) / vkarmn / zsqrt_CdN
	205	ztmp2 = zsqrt_Cd / zsqrt_CdN
	206
	207	ztmp1 = 0.5_wp + SIGN(0.5_wp,zeta_u) ! update stability
	208	zChN = 1.e-3_wp * zsqrt_CdN(18._wpztmp1 + 32.7_wp*(1._wp - ztmp1)) ! L&Y 2004 eq. (6c-6d)
	209	zCeN = 1.e-3_wp * (34.6_wp * zsqrt_CdN) ! L&Y 2004 eq. (6b)
	210
	211	Ch = MAX( zChNztmp2 / ( 1._wp + zChNztmp0 ) , Cx_min ) ! L&Y 2004 eq. (10b)
	212	Ce = MAX( zCeNztmp2 / ( 1._wp + zCeNztmp0 ) , Cx_min ) ! L&Y 2004 eq. (10c)
	213
[12377]	214	ENDIF
	215
[14072]	216	END DO !DO jit = 1, nbit
[12377]	217
[14072]	218	IF(PRESENT(CdN)) CdN(:,:) = zCdN(:,:)
	219	IF(PRESENT(CeN)) CeN(:,:) = zCeN(:,:)
	220	IF(PRESENT(ChN)) ChN(:,:) = zChN(:,:)
	221
[6723]	222	END SUBROUTINE turb_ncar
	223
	224
[14072]	225	FUNCTION cd_n10_ncar( pw10 )
	226	!!----------------------------------------------------------------------------------
[6723]	227	!! Estimate of the neutral drag coefficient at 10m as a function
	228	!! of neutral wind speed at 10m
	229	!!
[14072]	230	!! Origin: Large & Yeager 2008, Eq. (11)
[6723]	231	!!
[12377]	232	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	233	!!----------------------------------------------------------------------------------
	234	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s)
[14072]	235	REAL(wp), DIMENSION(jpi,jpj) :: cd_n10_ncar
[6723]	236	!
	237	INTEGER :: ji, jj ! dummy loop indices
	238	REAL(wp) :: zgt33, zw, zw6 ! local scalars
	239	!!----------------------------------------------------------------------------------
	240	!
[13460]	241	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[14072]	242	!
	243	zw = pw10(ji,jj)
	244	zw6 = zwzwzw
	245	zw6 = zw6*zw6
	246	!
	247	! When wind speed > 33 m/s => Cyclone conditions => special treatment
	248	zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1
	249	!
	250	cd_n10_ncar(ji,jj) = 1.e-3_wp * ( &
	251	& (1._wp - zgt33)( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wpzw6) & ! wind < 33 m/s
	252	& + zgt33 * 2.34_wp ) ! wind >= 33 m/s
	253	!
	254	cd_n10_ncar(ji,jj) = MAX( cd_n10_ncar(ji,jj), Cx_min )
	255	!
[12377]	256	END_2D
[6723]	257	!
[14072]	258	END FUNCTION cd_n10_ncar
[6723]	259
	260
[14072]	261	FUNCTION ch_n10_ncar( psqrtcdn10 , pstab )
[6723]	262	!!----------------------------------------------------------------------------------
[14072]	263	!! Estimate of the neutral heat transfer coefficient at 10m !!
	264	!! Origin: Large & Yeager 2008, Eq. (9) and (12)
	265
	266	!!----------------------------------------------------------------------------------
	267	REAL(wp), DIMENSION(jpi,jpj) :: ch_n10_ncar
	268	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
	269	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab ! stable ABL => 1 / unstable ABL => 0
	270	!!----------------------------------------------------------------------------------
	271	IF( ANY(pstab < -0.00001) .OR. ANY(pstab > 1.00001) ) THEN
	272	PRINT *, 'ERROR: ch_n10_ncar@mod_blk_ncar.f90: pstab ='
	273	PRINT *, pstab
	274	STOP
	275	END IF
	276	!
	277	ch_n10_ncar = MAX( 1.e-3_wp * psqrtcdn10( 18._wppstab + 32.7_wp*(1._wp - pstab) ) , Cx_min ) ! Eq. (9) & (12) Large & Yeager, 2008
	278	!
	279	END FUNCTION ch_n10_ncar
	280
	281	FUNCTION ce_n10_ncar( psqrtcdn10 )
	282	!!----------------------------------------------------------------------------------
	283	!! Estimate of the neutral heat transfer coefficient at 10m !!
	284	!! Origin: Large & Yeager 2008, Eq. (9) and (13)
	285	!!----------------------------------------------------------------------------------
	286	REAL(wp), DIMENSION(jpi,jpj) :: ce_n10_ncar
	287	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
	288	!!----------------------------------------------------------------------------------
	289	ce_n10_ncar = MAX( 1.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , Cx_min )
	290	!
	291	END FUNCTION ce_n10_ncar
	292
	293
	294	FUNCTION psi_m_ncar( pzeta )
	295	!!----------------------------------------------------------------------------------
[6723]	296	!! Universal profile stability function for momentum
[14072]	297	!! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
[12377]	298	!!
	299	!! pzeta : stability paramenter, z/L where z is altitude measurement
[6723]	300	!! and L is M-O length
	301	!!
[12377]	302	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	303	!!----------------------------------------------------------------------------------
[14072]	304	REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar
[12377]	305	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
[6723]	306	!
[12377]	307	INTEGER :: ji, jj ! dummy loop indices
[14072]	308	REAL(wp) :: zta, zx2, zx, zpsi_unst, zpsi_stab, zstab ! local scalars
[6723]	309	!!----------------------------------------------------------------------------------
[13460]	310	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[14072]	311	zta = pzeta(ji,jj)
	312	!
	313	zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5
	314	zx2 = MAX( zx2 , 1._wp )
	315	zx = SQRT(zx2) ! (1 - 16z)^0.25
	316	zpsi_unst = 2._wpLOG( (1._wp + zx )0.5_wp ) &
	317	& + LOG( (1._wp + zx2)*0.5_wp ) &
	318	& - 2._wpATAN(zx) + rpi0.5_wp
	319	!
	320	zpsi_stab = -5._wp*zta
	321	!
	322	zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
	323	!
	324	psi_m_ncar(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
	325	& + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
	326	!
	327	!
[12377]	328	END_2D
[14072]	329	END FUNCTION psi_m_ncar
[6723]	330
	331
[14072]	332	FUNCTION psi_h_ncar( pzeta )
[6723]	333	!!----------------------------------------------------------------------------------
	334	!! Universal profile stability function for temperature and humidity
[14072]	335	!! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
[6723]	336	!!
[12377]	337	!! pzeta : stability paramenter, z/L where z is altitude measurement
[6723]	338	!! and L is M-O length
	339	!!
[12377]	340	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	341	!!----------------------------------------------------------------------------------
[14072]	342	REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar
[6723]	343	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
	344	!
[12377]	345	INTEGER :: ji, jj ! dummy loop indices
[14072]	346	REAL(wp) :: zta, zx2, zpsi_unst, zpsi_stab, zstab ! local scalars
[6723]	347	!!----------------------------------------------------------------------------------
	348	!
[13460]	349	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[14072]	350	!
	351	zta = pzeta(ji,jj)
	352	!
	353	zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5
	354	zx2 = MAX( zx2 , 1._wp )
	355	zpsi_unst = 2._wpLOG( 0.5_wp(1._wp + zx2) )
	356	!
	357	zpsi_stab = -5._wp*zta
	358	!
	359	zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
	360	!
	361	psi_h_ncar(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
	362	& + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
	363	!
[12377]	364	END_2D
[14072]	365	END FUNCTION psi_h_ncar
[6723]	366
	367	!!======================================================================
	368	END MODULE sbcblk_algo_ncar

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 14205

Download in other formats: