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/2019/dev_ASINTER-01-05_merged/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 12015

Last change on this file since 12015 was 12015, checked in by gsamson, 4 years ago
dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK
Property svn:keywords set to `Id`
File size: 15.8 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
	7	!! * the effective bulk wind speed at 10m U_blk
	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
[12015]	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 oce ! ocean dynamics and tracers
	26	USE dom_oce ! ocean space and time domain
	27	USE phycst ! physical constants
[12015]	28	USE iom ! I/O manager library
	29	USE lib_mpp ! distribued memory computing library
	30	USE in_out_manager ! I/O manager
	31	USE prtctl ! Print control
[6723]	32	USE sbcwave, ONLY : cdn_wave ! wave module
[9570]	33	#if defined key_si3 \|\| defined key_cice
[6723]	34	USE sbc_ice ! Surface boundary condition: ice fields
	35	#endif
	36	USE lib_fortran ! to use key_nosignedzero
	37
[12015]	38	USE sbc_oce ! Surface boundary condition: ocean fields
	39	USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB
[6723]	40
	41	IMPLICIT NONE
	42	PRIVATE
	43
[12015]	44	PUBLIC :: TURB_NCAR ! called by sbcblk.F90
[6723]	45
[12015]	46	INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations
	47
[6723]	48	!!----------------------------------------------------------------------
	49	CONTAINS
	50
	51	SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
[9019]	52	& Cd, Ch, Ce, t_zu, q_zu, U_blk, &
	53	& Cdn, Chn, Cen )
[12015]	54	!!----------------------------------------------------------------------
[6723]	55	!! * ROUTINE turb_ncar *
	56	!!
	57	!! ** Purpose : Computes turbulent transfert coefficients of surface
	58	!! fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
	59	!! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
	60	!! Returns the effective bulk wind speed at 10m to be used in the bulk formulas
	61	!!
	62	!!
	63	!! INPUT :
	64	!! -------
	65	!! * zt : height for temperature and spec. hum. of air [m]
[12015]	66	!! * zu : height for wind speed (usually 10m) [m]
	67	!! * sst : bulk SST [K]
[6723]	68	!! * t_zt : potential air temperature at zt [K]
	69	!! * ssq : specific humidity at saturation at SST [kg/kg]
	70	!! * q_zt : specific humidity of air at zt [kg/kg]
[12015]	71	!! * U_zu : scalar wind speed at zu [m/s]
[6723]	72	!!
	73	!!
	74	!! OUTPUT :
	75	!! --------
	76	!! * Cd : drag coefficient
	77	!! * Ch : sensible heat coefficient
	78	!! * Ce : evaporation coefficient
	79	!! * t_zu : pot. air temperature adjusted at wind height zu [K]
	80	!! * q_zu : specific humidity of air // [kg/kg]
[12015]	81	!! * U_blk : bulk wind speed at zu [m/s]
	82	!!
	83	!!
	84	!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	85	!!----------------------------------------------------------------------------------
	86	REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m]
	87	REAL(wp), INTENT(in ) :: zu ! height for U_zu [m]
	88	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst ! sea surface temperature [Kelvin]
	89	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin]
	90	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq ! sea surface specific humidity [kg/kg]
[12015]	91	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg]
[6723]	92	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s]
	93	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau)
	94	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens)
	95	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat)
	96	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K]
	97	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg]
[12015]	98	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s]
[9019]	99	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients
[6723]	100	!
[12015]	101	INTEGER :: j_itt
	102	LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U
[6723]	103	!
[9125]	104	REAL(wp), DIMENSION(jpi,jpj) :: Cx_n10 ! 10m neutral latent/sensible coefficient
	105	REAL(wp), DIMENSION(jpi,jpj) :: sqrt_Cd_n10 ! root square of Cd_n10
	106	REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
	107	REAL(wp), DIMENSION(jpi,jpj) :: zpsi_h_u
	108	REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
	109	REAL(wp), DIMENSION(jpi,jpj) :: stab ! stability test integer
[6723]	110	!!----------------------------------------------------------------------------------
	111	!
	112	l_zt_equal_zu = .FALSE.
[12015]	113	IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision
[6723]	114
[12015]	115	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
[6723]	116
	117	!! First guess of stability:
[12015]	118	ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt
	119	stab = 0.5_wp + sign(0.5_wp,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable
[6723]	120
	121	!! Neutral coefficients at 10m:
	122	IF( ln_cdgw ) THEN ! wave drag case
[7753]	123	cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
	124	ztmp0 (:,:) = cdn_wave(:,:)
[6723]	125	ELSE
[12015]	126	ztmp0 = cd_neutral_10m( U_blk )
[6723]	127	ENDIF
	128
	129	sqrt_Cd_n10 = SQRT( ztmp0 )
	130
	131	!! Initializing transf. coeff. with their first guess neutral equivalents :
	132	Cd = ztmp0
[12015]	133	Ce = 1.e-3_wp( 34.6_wp sqrt_Cd_n10 )
	134	Ch = 1.e-3_wpsqrt_Cd_n10(18._wpstab + 32.7_wp(1._wp - stab))
[6723]	135	stab = sqrt_Cd_n10 ! Temporaty array !!! stab == SQRT(Cd)
[10190]	136
[12015]	137	IF( ln_cdgw ) THEN
	138	Cen = Ce
	139	Chn = Ch
	140	ENDIF
[6723]	141
	142	!! Initializing values at z_u with z_t values:
	143	t_zu = t_zt ; q_zu = q_zt
	144
[12015]	145	!! ITERATION BLOCK
	146	DO j_itt = 1, nb_itt
[6723]	147	!
	148	ztmp1 = t_zu - sst ! Updating air/sea differences
	149	ztmp2 = q_zu - ssq
	150
	151	! Updating turbulent scales : (L&Y 2004 eq. (7))
[12015]	152	ztmp0 = stabU_blk ! u (stab == SQRT(Cd))
	153	ztmp1 = Ch/stabztmp1 ! theta (stab == SQRT(Cd))
	154	ztmp2 = Ce/stabztmp2 ! q (stab == SQRT(Cd))
[6723]	155
	156	! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
[12015]	157	ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
	158
[6723]	159	!! Stability parameters :
[12015]	160	zeta_u = zu*ztmp0
	161	zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u )
[6723]	162	zpsi_h_u = psi_h( zeta_u )
	163
	164	!! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
	165	IF( .NOT. l_zt_equal_zu ) THEN
	166	!! Array 'stab' is free for the moment so using it to store 'zeta_t'
[12015]	167	stab = zt*ztmp0
	168	stab = SIGN( MIN(ABS(stab),10._wp), stab ) ! Temporaty array stab == zeta_t !!!
[6723]	169	stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again!
	170	t_zu = t_zt - ztmp1/vkarmnstab ! ztmp1 is still theta L&Y 2004 eq.(9b)
	171	q_zu = q_zt - ztmp2/vkarmnstab ! ztmp2 is still q L&Y 2004 eq.(9c)
[12015]	172	q_zu = max(0._wp, q_zu)
	173	ENDIF
[6723]	174
[12015]	175	! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
	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.
[6723]	179	ztmp2 = psi_m(zeta_u)
	180	IF( ln_cdgw ) THEN ! surface wave case
	181	stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 ) ! (stab == SQRT(Cd))
[10190]	182	Cd = stab * stab
[12015]	183	ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
[10190]	184	ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd))
[12015]	185	ztmp1 = 1._wp + Chn * ztmp0
[10190]	186	Ch = Chn * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b)
[12015]	187	ztmp1 = 1._wp + Cen * ztmp0
[10190]	188	Ce = Cen * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c)
	189
[6723]	190	ELSE
[12015]	191	! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
	192	! In very rare low-wind conditions, the old way of estimating the
	193	! neutral wind speed at 10m leads to a negative value that causes the code
	194	! to crash. To prevent this a threshold of 0.25m/s is imposed.
	195	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))
	196	ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10
	197	Cdn(:,:) = ztmp0
	198	sqrt_Cd_n10 = sqrt(ztmp0)
[6723]	199
[12015]	200	stab = 0.5_wp + sign(0.5_wp,zeta_u) ! update stability
	201	Cx_n10 = 1.e-3_wpsqrt_Cd_n10(18._wpstab + 32.7_wp(1._wp - stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10)
	202	Chn(:,:) = Cx_n10
[6723]	203
[12015]	204	!! Update of transfer coefficients:
	205	ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
	206	Cd = ztmp0 / ( ztmp1*ztmp1 )
	207	stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
[6723]	208
[12015]	209	ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
	210	ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd))
	211	ztmp1 = 1._wp + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10)
	212	Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b)
[6723]	213
[12015]	214	Cx_n10 = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10
	215	Cen(:,:) = Cx_n10
	216	ztmp1 = 1._wp + Cx_n10*ztmp0
	217	Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c)
	218	ENDIF
	219
	220	END DO !DO j_itt = 1, nb_itt
	221
[6723]	222	END SUBROUTINE turb_ncar
	223
	224
	225	FUNCTION cd_neutral_10m( pw10 )
	226	!!----------------------------------------------------------------------------------
	227	!! Estimate of the neutral drag coefficient at 10m as a function
	228	!! of neutral wind speed at 10m
	229	!!
	230	!! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
	231	!!
[12015]	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)
	235	REAL(wp), DIMENSION(jpi,jpj) :: cd_neutral_10m
	236	!
	237	INTEGER :: ji, jj ! dummy loop indices
	238	REAL(wp) :: zgt33, zw, zw6 ! local scalars
	239	!!----------------------------------------------------------------------------------
	240	!
	241	DO jj = 1, jpj
	242	DO ji = 1, jpi
	243	!
	244	zw = pw10(ji,jj)
	245	zw6 = zwzwzw
	246	zw6 = zw6*zw6
	247	!
	248	! When wind speed > 33 m/s => Cyclone conditions => special treatment
[12015]	249	zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1
[6723]	250	!
[12015]	251	cd_neutral_10m(ji,jj) = 1.e-3_wp * ( &
	252	& (1._wp - zgt33)( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wpzw6) & ! wind < 33 m/s
	253	& + zgt33 * 2.34_wp ) ! wind >= 33 m/s
[6723]	254	!
[12015]	255	cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp)
[6723]	256	!
	257	END DO
	258	END DO
	259	!
	260	END FUNCTION cd_neutral_10m
	261
	262
	263	FUNCTION psi_m( pzeta )
	264	!!----------------------------------------------------------------------------------
	265	!! Universal profile stability function for momentum
	266	!! !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
[12015]	267	!!
	268	!! pzeta : stability paramenter, z/L where z is altitude measurement
[6723]	269	!! and L is M-O length
	270	!!
[12015]	271	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	272	!!----------------------------------------------------------------------------------
[12015]	273	REAL(wp), DIMENSION(jpi,jpj) :: psi_m
	274	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
[6723]	275	!
[12015]	276	INTEGER :: ji, jj ! dummy loop indices
[6723]	277	REAL(wp) :: zx2, zx, zstab ! local scalars
	278	!!----------------------------------------------------------------------------------
	279	DO jj = 1, jpj
	280	DO ji = 1, jpi
[12015]	281	zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
	282	zx2 = MAX( zx2 , 1._wp )
[6723]	283	zx = SQRT( zx2 )
[12015]	284	zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
[6723]	285	!
[12015]	286	psi_m(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable
	287	& + (1._wp - zstab) * (2._wpLOG((1._wp + zx)0.5_wp) & ! Unstable
	288	& + LOG((1._wp + zx2)0.5_wp) - 2._wpATAN(zx) + rpi*0.5_wp) ! "
[6723]	289	!
	290	END DO
	291	END DO
	292	END FUNCTION psi_m
	293
	294
	295	FUNCTION psi_h( pzeta )
	296	!!----------------------------------------------------------------------------------
	297	!! Universal profile stability function for temperature and humidity
	298	!! !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
	299	!!
[12015]	300	!! pzeta : stability paramenter, z/L where z is altitude measurement
[6723]	301	!! and L is M-O length
	302	!!
[12015]	303	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	304	!!----------------------------------------------------------------------------------
[12015]	305	REAL(wp), DIMENSION(jpi,jpj) :: psi_h
[6723]	306	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
	307	!
[12015]	308	INTEGER :: ji, jj ! dummy loop indices
[6723]	309	REAL(wp) :: zx2, zstab ! local scalars
	310	!!----------------------------------------------------------------------------------
	311	!
	312	DO jj = 1, jpj
	313	DO ji = 1, jpi
[12015]	314	zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
	315	zx2 = MAX( zx2 , 1._wp )
	316	zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
[6723]	317	!
[12015]	318	psi_h(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable
	319	& + (1._wp - zstab) * (2._wpLOG( (1._wp + zx2)0.5_wp )) ! Unstable
[6723]	320	!
	321	END DO
	322	END DO
	323	END FUNCTION psi_h
	324
	325	!!======================================================================
	326	END MODULE sbcblk_algo_ncar

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

Download in other formats: