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/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 13673

Last change on this file since 13673 was 13017, checked in by laurent, 4 years ago
Bugfix to prevent "TURB_NCAR()" to return inconsistsent Cd and CdN when "ln_wave==T".
Property svn:keywords set to `Id`
File size: 21.0 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
[13017]	13	!! (https://github.com/brodeau/aerobulk/)
[6723]	14	!!
[13017]	15	!! L. Brodeau, 2020
[6723]	16	!!=====================================================================
[13017]	17	!! History : 4.0 ! 2020-06 (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
	28	USE sbc_oce ! Surface boundary condition: ocean fields
	29	USE sbcwave, ONLY : cdn_wave ! wave module
[9570]	30	#if defined key_si3 \|\| defined key_cice
[6723]	31	USE sbc_ice ! Surface boundary condition: ice fields
	32	#endif
	33	!
	34	USE iom ! I/O manager library
	35	USE lib_mpp ! distribued memory computing library
	36	USE in_out_manager ! I/O manager
	37	USE prtctl ! Print control
	38	USE lib_fortran ! to use key_nosignedzero
	39
	40
	41	IMPLICIT NONE
	42	PRIVATE
	43
[13017]	44	PUBLIC :: TURB_NCAR ! called by sbcblk.F90
[6723]	45
[6727]	46	! ! NCAR own values for given constants:
[6723]	47	REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature...
[13017]	48
	49	INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations
	50
[6723]	51	!!----------------------------------------------------------------------
	52	CONTAINS
	53
	54	SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
[13017]	55	& Cd, Ch, Ce, t_zu, q_zu, Ub, &
	56	& CdN, ChN, CeN )
	57	!!----------------------------------------------------------------------
[6723]	58	!! * ROUTINE turb_ncar *
	59	!!
	60	!! ** Purpose : Computes turbulent transfert coefficients of surface
	61	!! fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
	62	!! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
[13017]	63	!! Returns the effective bulk wind speed at zu to be used in the bulk formulas
[6723]	64	!!
	65	!! INPUT :
	66	!! -------
	67	!! * zt : height for temperature and spec. hum. of air [m]
[13017]	68	!! * zu : height for wind speed (usually 10m) [m]
	69	!! * sst : bulk SST [K]
[6723]	70	!! * t_zt : potential air temperature at zt [K]
	71	!! * ssq : specific humidity at saturation at SST [kg/kg]
	72	!! * q_zt : specific humidity of air at zt [kg/kg]
[13017]	73	!! * U_zu : scalar wind speed at zu [m/s]
[6723]	74	!!
	75	!! OUTPUT :
	76	!! --------
	77	!! * Cd : drag coefficient
	78	!! * Ch : sensible heat coefficient
	79	!! * Ce : evaporation coefficient
	80	!! * t_zu : pot. air temperature adjusted at wind height zu [K]
	81	!! * q_zu : specific humidity of air // [kg/kg]
[13017]	82	!! * Ub : bulk wind speed at zu [m/s]
	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]
[13017]	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]
[13017]	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
[6723]	100	!
[13017]	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
[13017]	105	REAL(wp), DIMENSION(jpi,jpj) :: sqrtCdN10 ! root square of Cd_n10
[9125]	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
[13017]	109	REAL(wp), DIMENSION(jpi,jpj) :: sqrtCd
[6723]	110	!!----------------------------------------------------------------------------------
	111
[13017]	112	l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp )
[6723]	113
[13017]	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
[6723]	115
[13017]	116	!! Neutral drag coefficient at zu:
[6723]	117	IF( ln_cdgw ) THEN ! wave drag case
[13017]	118	CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp )
[6723]	119	ELSE
[13017]	120	CdN = CD_N10_NCAR( Ub )
[6723]	121	ENDIF
[13017]	122	sqrtCdN10 = SQRT( CdN )
[6723]	123
	124	!! Initializing transf. coeff. with their first guess neutral equivalents :
[13017]	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
[6723]	130
	131	!! Initializing values at z_u with z_t values:
[13017]	132	t_zu = t_zt
	133	q_zu = q_zt
[6723]	134
[13017]	135	!! ITERATION BLOCK
	136	DO j_itt = 1, nb_itt
[6723]	137	!
	138	ztmp1 = t_zu - sst ! Updating air/sea differences
	139	ztmp2 = q_zu - ssq
	140
[13017]	141	! Updating turbulent scales : (L&Y 2004 Eq. (7))
	142	ztmp0 = sqrtCdUb ! u
	143	ztmp1 = Ch/sqrtCdztmp1 ! theta
	144	ztmp2 = Ce/sqrtCdztmp2 ! q
[6723]	145
[13017]	146	! Estimate the inverse of Obukov length (1/L) at height zu:
	147	ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
[6723]	148
	149	!! Stability parameters :
[13017]	150	zeta_u = zu*ztmp0
	151	zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u )
[6723]	152
[13017]	153	!! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b-9c))
[6723]	154	IF( .NOT. l_zt_equal_zu ) THEN
[13017]	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/vkarmnztmp0 ! ztmp1 is still theta L&Y 2004 Eq. (9b)
	159	!!
	160	q_zu = q_zt - ztmp2/vkarmnztmp0 ! ztmp2 is still q L&Y 2004 Eq. (9c)
	161	q_zu = MAX(0._wp, q_zu)
[6723]	162	END IF
	163
[13017]	164	! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)...
	165	! In very rare low-wind 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))
[10190]	170
[13017]	171	IF( ln_cdgw ) THEN ! wave drag case
	172	CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp )
[6723]	173	ELSE
[13017]	174	CdN = CD_N10_NCAR(ztmp0) ! Cd_n10
	175	END IF
	176	sqrtCdN10 = SQRT(CdN)
[6723]	177
[13017]	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.1E-3_wp )
	181	sqrtCd = SQRT( Cd )
[6723]	182
[13017]	183	ztmp0 = ( LOG(zu/10._wp) - psi_h_ncar(zeta_u) ) / vkarmn / sqrtCdN10
	184	ztmp2 = sqrtCd / sqrtCdN10
[6723]	185
[13017]	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.1E-3_wp ) ! L&Y 2004 Eq. (10b)
[6723]	190
[13017]	191	CeN = CE_N10_NCAR( sqrtCdN10 )
	192	ztmp1 = 1._wp + CeN*ztmp0
	193	Ce = MAX( CeN*ztmp2 / ztmp1 , 0.1E-3_wp ) ! L&Y 2004 Eq. (10c)
	194
	195	END DO !DO j_itt = 1, nb_itt
	196
[6723]	197	END SUBROUTINE turb_ncar
	198
	199
[13017]	200	FUNCTION CD_N10_NCAR( pw10 )
	201	!!----------------------------------------------------------------------------------
[6723]	202	!! Estimate of the neutral drag coefficient at 10m as a function
	203	!! of neutral wind speed at 10m
	204	!!
[13017]	205	!! Origin: Large & Yeager 2008, Eq. (11)
[6723]	206	!!
[13017]	207	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	208	!!----------------------------------------------------------------------------------
	209	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s)
[13017]	210	REAL(wp), DIMENSION(jpi,jpj) :: CD_N10_NCAR
[6723]	211	!
	212	INTEGER :: ji, jj ! dummy loop indices
	213	REAL(wp) :: zgt33, zw, zw6 ! local scalars
	214	!!----------------------------------------------------------------------------------
	215	!
	216	DO jj = 1, jpj
	217	DO ji = 1, jpi
	218	!
	219	zw = pw10(ji,jj)
	220	zw6 = zwzwzw
	221	zw6 = zw6*zw6
	222	!
	223	! When wind speed > 33 m/s => Cyclone conditions => special treatment
[13017]	224	zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1
[6723]	225	!
[13017]	226	CD_N10_NCAR(ji,jj) = 1.e-3_wp * ( &
	227	& (1._wp - zgt33)( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wpzw6) & ! wind < 33 m/s
	228	& + zgt33 * 2.34_wp ) ! wind >= 33 m/s
[6723]	229	!
[13017]	230	CD_N10_NCAR(ji,jj) = MAX( CD_N10_NCAR(ji,jj), 0.1E-3_wp )
[6723]	231	!
	232	END DO
	233	END DO
	234	!
[13017]	235	END FUNCTION CD_N10_NCAR
[6723]	236
	237
[13017]	238
	239	FUNCTION CH_N10_NCAR( psqrtcdn10 , pstab )
[6723]	240	!!----------------------------------------------------------------------------------
[13017]	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.e-3_wp * psqrtcdn10( 18._wppstab + 32.7_wp*(1._wp - pstab) ) , 0.1E-3_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.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , 0.1E-3_wp )
	263	!
	264	END FUNCTION CE_N10_NCAR
	265
	266
	267	FUNCTION psi_m_ncar( pzeta )
	268	!!----------------------------------------------------------------------------------
[6723]	269	!! Universal profile stability function for momentum
[13017]	270	!! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
	271	!!
	272	!! pzeta : stability paramenter, z/L where z is altitude measurement
[6723]	273	!! and L is M-O length
	274	!!
[13017]	275	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	276	!!----------------------------------------------------------------------------------
[13017]	277	REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar
	278	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
[6723]	279	!
[13017]	280	INTEGER :: ji, jj ! dummy loop indices
	281	REAL(wp) :: zzeta, zx2, zx, zpsi_unst, zpsi_stab, zstab ! local scalars
[6723]	282	!!----------------------------------------------------------------------------------
	283	DO jj = 1, jpj
	284	DO ji = 1, jpi
[13017]	285
	286	zzeta = pzeta(ji,jj)
[6723]	287	!
[13017]	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._wpLOG( (1._wp + zx )0.5_wp ) &
	292	& + LOG( (1._wp + zx2)*0.5_wp ) &
	293	& - 2._wpATAN(zx) + rpi0.5_wp
[6723]	294	!
[13017]	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
	301	!
[6723]	302	END DO
	303	END DO
[13017]	304	END FUNCTION psi_m_ncar
[6723]	305
	306
[13017]	307	FUNCTION psi_h_ncar( pzeta )
[6723]	308	!!----------------------------------------------------------------------------------
	309	!! Universal profile stability function for temperature and humidity
[13017]	310	!! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
[6723]	311	!!
[13017]	312	!! pzeta : stability paramenter, z/L where z is altitude measurement
[6723]	313	!! and L is M-O length
	314	!!
[13017]	315	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	316	!!----------------------------------------------------------------------------------
[13017]	317	REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar
[6723]	318	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
	319	!
[13017]	320	INTEGER :: ji, jj ! dummy loop indices
	321	REAL(wp) :: zzeta, zx2, zpsi_unst, zpsi_stab, zstab ! local scalars
[6723]	322	!!----------------------------------------------------------------------------------
	323	!
	324	DO jj = 1, jpj
	325	DO ji = 1, jpi
	326	!
[13017]	327	zzeta = pzeta(ji,jj)
[6723]	328	!
[13017]	329	zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) ) ! (1 -16z)^0.5
	330	zx2 = MAX( zx2 , 1._wp )
	331	zpsi_unst = 2._wpLOG( 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
	339	!
[6723]	340	END DO
	341	END DO
[13017]	342	END FUNCTION psi_h_ncar
	343
	344
	345
	346
	347	FUNCTION UN10_from_CD( pzu, pUb, pCd, ppsi )
	348	!!----------------------------------------------------------------------------------
	349	!! Provides the neutral-stability 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.
[6723]	379	!
[13017]	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) = gravvkarmn( pts(ji,jj)zqa + rctv0ptha(ji,jj)*pqs(ji,jj) ) &
	395	& / MAX( pus(ji,jj)pus(ji,jj)ptha(ji,jj)*zqa , 1.E-9_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
[6723]	403
[13017]	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 non-neutral) 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 NEUTRAL-STABILITY 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 neutral-stability CdN) :
	416	z0_from_Cd = pzu * EXP( - ( vkarmn/SQRT(pCd(:,:)) + ppsi(:,:) ) ) !LB: ok, double-checked!
	417	ELSE
	418	!! Cd provided is the neutral-stability Cd, aka CdN :
	419	z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) ) !LB: ok, double-checked!
	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
	429
[6723]	430	!!======================================================================
	431	END MODULE sbcblk_algo_ncar

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

Download in other formats: