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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcblk_algo_ncar.F90 in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 13655

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

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

Download in other formats: