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_coare.F90 in NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90 @ 11266

Last change on this file since 11266 was 11266, checked in by laurent, 5 years ago
LB: CSWL param now uses NEMO time step (rdt) and previous-t value of t_skin as first guess.
Property svn:keywords set to `Id`
File size: 21.9 KB

Rev	Line
[6723]	1	MODULE sbcblk_algo_coare
	2	!!======================================================================
[11111]	3	!! * MODULE sbcblk_algo_coare *
	4	!! Computes:
[6723]	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	!!
[11111]	10	!! Using the bulk formulation/param. of COARE v3, Fairall et al., 2003
[6723]	11	!!
	12	!! Routine turb_coare maintained and developed in AeroBulk
[11111]	13	!! (https://github.com/brodeau/aerobulk)
[6723]	14	!!
[11215]	15	!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
[6727]	16	!!----------------------------------------------------------------------
[11111]	17	!! History : 4.0 ! 2016-02 (L.Brodeau) Original code
	18	!!----------------------------------------------------------------------
[6723]	19
	20	!!----------------------------------------------------------------------
	21	!! turb_coare : 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	!!----------------------------------------------------------------------
[11111]	25	USE oce ! ocean dynamics and tracers
	26	USE dom_oce ! ocean space and time domain
	27	USE phycst ! physical constants
[11215]	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
[11111]	32	USE sbcwave, ONLY : cdn_wave ! wave module
[9570]	33	#if defined key_si3 \|\| defined key_cice
[11111]	34	USE sbc_ice ! Surface boundary condition: ice fields
[6723]	35	#endif
[11215]	36	USE lib_fortran ! to use key_nosignedzero
[6723]	37
[11215]	38	USE sbc_oce ! Surface boundary condition: ocean fields
	39	USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB
	40	USE sbcblk_skin ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
	41
[6723]	42	IMPLICIT NONE
	43	PRIVATE
	44
[6727]	45	PUBLIC :: TURB_COARE ! called by sbcblk.F90
[6723]	46
[9019]	47	! !! COARE own values for given constants:
	48	REAL(wp), PARAMETER :: zi0 = 600._wp ! scale height of the atmospheric boundary layer...
	49	REAL(wp), PARAMETER :: Beta0 = 1.250_wp ! gustiness parameter
[6723]	50
[11111]	51	INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations
	52
[9019]	53	!!----------------------------------------------------------------------
[6723]	54	CONTAINS
	55
[11215]	56	SUBROUTINE turb_coare( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
[9019]	57	& Cd, Ch, Ce, t_zu, q_zu, U_blk, &
[11215]	58	& Cdn, Chn, Cen, &
[11266]	59	& Qsw, rad_lw, slp, Tsk_b )
[6723]	60	!!----------------------------------------------------------------------
	61	!! * ROUTINE turb_coare *
	62	!!
	63	!! ** Purpose : Computes turbulent transfert coefficients of surface
	64	!! fluxes according to Fairall et al. (2003)
	65	!! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
[11111]	66	!! Returns the effective bulk wind speed at 10m to be used in the bulk formulas
[6723]	67	!!
[11215]	68	!! Applies the cool-skin warm-layer correction of the SST to T_s
	69	!! if the net shortwave flux at the surface (Qsw), the downwelling longwave
	70	!! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
	71	!! are provided as (optional) arguments!
[6723]	72	!!
	73	!! INPUT :
	74	!! -------
	75	!! * zt : height for temperature and spec. hum. of air [m]
	76	!! * zu : height for wind speed (generally 10m) [m]
	77	!! * U_zu : scalar wind speed at 10m [m/s]
	78	!! * t_zt : potential air temperature at zt [K]
	79	!! * q_zt : specific humidity of air at zt [kg/kg]
	80	!!
[11215]	81	!! INPUT/OUTPUT:
	82	!! -------------
[11266]	83	!! * T_s : always "bulk SST" as input [K]
	84	!! -> unchanged "bulk SST" as output if CSWL not used [K]
	85	!! -> skin temperature as output if CSWL used [K]
	86	!!
[11215]	87	!! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg]
	88	!! -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True)
	89	!! -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False)
[6723]	90	!!
[11215]	91	!! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!):
	92	!! ---------------
	93	!! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2]
	94	!! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2]
	95	!! * slp : sea-level pressure [Pa]
[11266]	96	!! * Tsk_b : estimate of skin temperature at previous time-step [K]
[11215]	97	!!
[6723]	98	!! OUTPUT :
	99	!! --------
	100	!! * Cd : drag coefficient
	101	!! * Ch : sensible heat coefficient
	102	!! * Ce : evaporation coefficient
	103	!! * t_zu : pot. air temperature adjusted at wind height zu [K]
	104	!! * q_zu : specific humidity of air // [kg/kg]
[11111]	105	!! * U_blk : bulk wind speed at 10m [m/s]
	106	!!
	107	!!
[11215]	108	!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
[11111]	109	!!----------------------------------------------------------------------------------
[6727]	110	REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m]
	111	REAL(wp), INTENT(in ) :: zu ! height for U_zu [m]
[11215]	112	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin]
[6727]	113	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin]
[11215]	114	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg]
[6727]	115	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg]
	116	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s]
[6723]	117	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau)
	118	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens)
	119	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat)
[6727]	120	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K]
	121	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg]
[6723]	122	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind at 10m [m/s]
[9019]	123	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients
[6727]	124	!
[11215]	125	REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2]
	126	REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2]
	127	REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa]
[11266]	128	REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Tsk_b ! [Pa]
[11215]	129	!
[6723]	130	INTEGER :: j_itt
[11111]	131	LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U
[11182]	132	!
[9125]	133	REAL(wp), DIMENSION(jpi,jpj) :: &
[6723]	134	& u_star, t_star, q_star, &
	135	& dt_zu, dq_zu, &
	136	& znu_a, & !: Nu_air, Viscosity of air
	137	& z0, z0t
[9125]	138	REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
	139	REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
	140	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t ! stability parameter at height zt
[6727]	141	!
[11215]	142	! Cool skin:
	143	LOGICAL :: l_use_skin = .FALSE.
	144	REAL(wp), DIMENSION(jpi,jpj) :: zsst ! to back up the initial bulk SST
	145	!!----------------------------------------------------------------------------------
	146	!
	147	! Cool skin ?
	148	IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
	149	l_use_skin = .TRUE.
	150	END IF
	151	IF (lwp) PRINT , ' ** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
	152
[6723]	153	l_zt_equal_zu = .FALSE.
[11215]	154	IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision
[6723]	155
[9125]	156	IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) )
[6723]	157
[11215]	158	!! Initialization for cool skin:
[11266]	159	zsst = T_s ! save the bulk SST
[11215]	160	IF( l_use_skin ) THEN
[11266]	161	! First guess for skin temperature:
	162	IF( PRESENT(Tsk_b) ) THEN
	163	T_s = Tsk_b
	164	ELSE
	165	T_s = T_s - 0.25 ! sst - 0.25
	166	END IF
[11215]	167	q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
	168	END IF
	169
[6723]	170	!! First guess of temperature and humidity at height zu:
[11215]	171	t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions...
[11111]	172	q_zu = MAX( q_zt , 1.e-6_wp ) ! "
[6723]	173
	174	!! Pot. temp. difference (and we don't want it to be 0!)
[11215]	175	dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
	176	dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
[6723]	177
[11111]	178	znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
[6723]	179
[11215]	180	ztmp2 = 0.5_wp*0.5_wp ! initial guess for wind gustiness contribution
[6723]	181	U_blk = SQRT(U_zu*U_zu + ztmp2)
	182
[11215]	183	ztmp2 = 10000._wp ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
[6723]	184	ztmp0 = LOG(zu*ztmp2)
	185	ztmp1 = LOG(10.*ztmp2)
[11215]	186	u_star = 0.035_wpU_blkztmp1/ztmp0 ! (u* = 0.035*Un10)
[6723]	187
[11215]	188	z0 = alfa_charn(U_blk)u_staru_star/grav + 0.11_wp*znu_a/u_star
	189	z0 = MIN(ABS(z0), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
	190	z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
	191	z0t = MIN(ABS(z0t), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
[6723]	192
	193	ztmp2 = vkarmn/ztmp0
	194	Cd = ztmp2*ztmp2 ! first guess of Cd
	195
[11111]	196	ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
[6723]	197
[11215]	198	ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
[11111]	199
[11209]	200	!! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
[11111]	201	ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
[6723]	202	ztmp0 = ztmp0*ztmp2
[11215]	203	zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi00.004_wpBeta0**3)))) & ! BRN < 0
	204	& + ztmp1 * (ztmp0(1._wp + 27._wp/9._wpztmp2/ztmp0)) ! BRN > 0
	205	!#LB: should make sure that the "ztmp0" of "27./9.ztmp2/ztmp0" is "ztmp0ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
[6723]	206
	207	!! First guess M-O stability dependent scaling params.(u,t,q*) to estimate z0 and z/L
[11111]	208	ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
[6723]	209
	210	u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
	211	t_star = dt_zu*ztmp0
	212	q_star = dq_zu*ztmp0
	213
[11266]	214	! What needs to be done if zt /= zu:
[6723]	215	IF( .NOT. l_zt_equal_zu ) THEN
[11111]	216	!! First update of values at zu (or zt for wind)
[6723]	217	zeta_t = zt*zeta_u/zu
	218	ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
[11111]	219	ztmp1 = LOG(zt/zu) + ztmp0
[6723]	220	t_zu = t_zt - t_star/vkarmn*ztmp1
	221	q_zu = q_zt - q_star/vkarmn*ztmp1
[11215]	222	q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
[11111]	223	!
[11215]	224	dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
	225	dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
[6723]	226	END IF
	227
	228	!! ITERATION BLOCK
	229	DO j_itt = 1, nb_itt
	230
	231	!!Inverse of Monin-Obukov length (1/L) :
	232	ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star) ! 1/L == 1/[Monin-Obukhov length]
[11111]	233	ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO
[6723]	234
	235	ztmp1 = u_staru_star ! u^2
	236
	237	!! Update wind at 10m taking into acount convection-related wind gustiness:
	238	! Ug = Betaw (Beta = 1.25, Fairall et al. 2003, Eq.8):
[11111]	239	ztmp2 = Beta0Beta0ztmp1(MAX(-zi0ztmp0/vkarmn,0._wp))**(2./3.) ! => ztmp2 == Ug^2
[6723]	240	!! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
[11111]	241	U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed
[6723]	242	! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
	243
	244	!! Updating Charnock parameter, increases with the wind (Fairall et al., 2003 p. 577-578)
	245	ztmp2 = alfa_charn(U_blk) ! alpha Charnock parameter
	246
	247	!! Roughness lengthes z0, z0t (z0q = z0t) :
[11215]	248	z0 = ztmp2ztmp1/grav + 0.11_wpznu_a/u_star ! Roughness length (eq.6)
	249	ztmp1 = z0*u_star/znu_a ! Re_r: roughness Reynolds number
[11111]	250	z0t = MIN( 1.1E-4_wp , 5.5E-5_wpztmp1*(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28)
[6723]	251
	252	!! Stability parameters:
[11111]	253	zeta_u = zu*ztmp0
	254	zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
[6723]	255	IF( .NOT. l_zt_equal_zu ) THEN
[11111]	256	zeta_t = zt*ztmp0
	257	zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
[6723]	258	END IF
	259
	260	!! Turbulent scales at zu=10m :
	261	ztmp0 = psi_h_coare(zeta_u)
[11111]	262	ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
[6723]	263
	264	t_star = dt_zu*ztmp1
	265	q_star = dq_zu*ztmp1
[11111]	266	u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
[6723]	267
	268	IF( .NOT. l_zt_equal_zu ) THEN
	269	! What's need to be done if zt /= zu
	270	!! Re-updating temperature and humidity at zu :
	271	ztmp2 = ztmp0 - psi_h_coare(zeta_t)
	272	ztmp1 = log(zt/zu) + ztmp2
	273	t_zu = t_zt - t_star/vkarmn*ztmp1
	274	q_zu = q_zt - q_star/vkarmn*ztmp1
	275	END IF
	276
[11215]	277	!! SKIN related part
	278	!! -----------------
	279	IF( l_use_skin ) THEN
	280	!! compute transfer coefficients at zu : lolo: verifier...
	281	ztmp0 = u_star/U_blk
	282	Ch = ztmp0*t_star/dt_zu
	283	Ce = ztmp0*q_star/dq_zu
	284	! Non-Solar heat flux to the ocean:
	285	ztmp1 = U_blkMAX(rho_air(t_zu, q_zu, slp), 1._wp) ! rhoU10
	286	ztmp2 = T_s*T_s
[11266]	287	ztmp1 = ztmp1 * ( CeL_vap(T_s)(q_zu - q_s) + Chcp_air(q_zu)(t_zu - T_s) ) & ! Total turb. heat flux
	288	& + rad_lw - emiss_wstefanztmp2*ztmp2 ! Net longwave flux
	289	!! => "ztmp1" is the net non-solar surface heat flux !
[11215]	290	!! Updating the values of the skin temperature T_s and q_s :
	291	CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
	292	q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! 200 -> just to avoid numerics problem on masked regions if silly values are given
	293	END IF
	294
	295	IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
	296	dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
	297	dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
	298	END IF
	299
	300	END DO !DO j_itt = 1, nb_itt
	301
[6723]	302	! compute transfer coefficients at zu :
	303	ztmp0 = u_star/U_blk
	304	Cd = ztmp0*ztmp0
	305	Ch = ztmp0*t_star/dt_zu
	306	Ce = ztmp0*q_star/dq_zu
	307	!
[9019]	308	ztmp1 = zu + z0
	309	Cdn = vkarmnvkarmn / (log(ztmp1/z0 )log(ztmp1/z0 ))
	310	Chn = vkarmnvkarmn / (log(ztmp1/z0t)log(ztmp1/z0t))
	311	Cen = Chn
	312	!
[9125]	313	IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
[9124]	314	!
[6723]	315	END SUBROUTINE turb_coare
	316
	317
	318	FUNCTION alfa_charn( pwnd )
	319	!!-------------------------------------------------------------------
	320	!! Compute the Charnock parameter as a function of the wind speed
	321	!!
	322	!! (Fairall et al., 2003 p.577-578)
	323	!!
	324	!! Wind below 10 m/s : alfa = 0.011
	325	!! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
	326	!! Wind greater than 18 m/s : alfa = 0.018
	327	!!
[11215]	328	!! Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	329	!!-------------------------------------------------------------------
	330	REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
	331	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed
	332	!
	333	INTEGER :: ji, jj ! dummy loop indices
	334	REAL(wp) :: zw, zgt10, zgt18
[6727]	335	!!-------------------------------------------------------------------
[6723]	336	!
	337	DO jj = 1, jpj
	338	DO ji = 1, jpi
	339	!
	340	zw = pwnd(ji,jj) ! wind speed
	341	!
	342	! Charnock's constant, increases with the wind :
[11111]	343	zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1
	344	zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
[6723]	345	!
	346	alfa_charn(ji,jj) = (1. - zgt10)*0.011 & ! wind is lower than 10 m/s
	347	& + zgt10((1. - zgt18)(0.011 + (0.018 - 0.011) &
	348	& (zw - 10.)/(18. - 10.)) + zgt18( 0.018 ) ) ! Hare et al. (1999)
	349	!
	350	END DO
	351	END DO
	352	!
	353	END FUNCTION alfa_charn
	354
	355	FUNCTION psi_m_coare( pzeta )
	356	!!----------------------------------------------------------------------------------
	357	!! ** Purpose: compute the universal profile stability function for momentum
	358	!! COARE 3.0, Fairall et al. 2003
	359	!! pzeta : stability paramenter, z/L where z is altitude
	360	!! measurement and L is M-O length
	361	!! Stability function for wind speed and scalars matching Kansas and free
	362	!! convection forms with weighting f convective form, follows Fairall et
	363	!! al (1996) with profile constants from Grachev et al (2000) BLM stable
	364	!! form from Beljaars and Holtslag (1991)
	365	!!
[11215]	366	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	367	!!----------------------------------------------------------------------------------
	368	REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
	369	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
	370	!
	371	INTEGER :: ji, jj ! dummy loop indices
	372	REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
	373	!!----------------------------------------------------------------------------------
	374	!
	375	DO jj = 1, jpj
	376	DO ji = 1, jpi
	377	!
	378	zta = pzeta(ji,jj)
	379	!
	380	zphi_m = ABS(1. - 15.zta)*.25 !!Kansas unstable
	381	!
	382	zpsi_k = 2.LOG((1. + zphi_m)/2.) + LOG((1. + zphi_mzphi_m)/2.) &
	383	& - 2.ATAN(zphi_m) + 0.5rpi
	384	!
	385	zphi_c = ABS(1. - 10.15zta)*.3333 !!Convective
	386	!
	387	zpsi_c = 1.5LOG((1. + zphi_c + zphi_czphi_c)/3.) &
	388	& - 1.7320508ATAN((1. + 2.zphi_c)/1.7320508) + 1.813799447
	389	!
	390	zf = zta*zta
	391	zf = zf/(1. + zf)
[11111]	392	zc = MIN(50._wp, 0.35_wp*zta)
	393	zstab = 0.5 + SIGN(0.5_wp, zta)
[6723]	394	!
	395	psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)zpsi_k + zfzpsi_c ) & ! (zta < 0)
	396	& - zstab * ( 1. + 1.*zta & ! (zta > 0)
	397	& + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 ) ! "
	398	!
	399	END DO
	400	END DO
	401	!
	402	END FUNCTION psi_m_coare
	403
	404
	405	FUNCTION psi_h_coare( pzeta )
	406	!!---------------------------------------------------------------------
	407	!! Universal profile stability function for temperature and humidity
	408	!! COARE 3.0, Fairall et al. 2003
	409	!!
	410	!! pzeta : stability paramenter, z/L where z is altitude measurement
	411	!! and L is M-O length
	412	!!
	413	!! Stability function for wind speed and scalars matching Kansas and free
	414	!! convection forms with weighting f convective form, follows Fairall et
	415	!! al (1996) with profile constants from Grachev et al (2000) BLM stable
	416	!! form from Beljaars and Holtslag (1991)
	417	!!
[11215]	418	!! Author: L. Brodeau, June 2016 / AeroBulk
[11111]	419	!! (https://github.com/brodeau/aerobulk/)
[6723]	420	!!----------------------------------------------------------------
	421	!!
	422	REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
	423	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
	424	!
	425	INTEGER :: ji, jj ! dummy loop indices
	426	REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
	427	!
	428	DO jj = 1, jpj
	429	DO ji = 1, jpi
	430	!
	431	zta = pzeta(ji,jj)
	432	!
	433	zphi_h = (ABS(1. - 15.zta)).5 !! Kansas unstable (zphi_h = zphi_m*2 when unstable, zphi_m when stable)
	434	!
	435	zpsi_k = 2.*LOG((1. + zphi_h)/2.)
	436	!
	437	zphi_c = (ABS(1. - 34.15zta))*.3333 !! Convective
	438	!
	439	zpsi_c = 1.5LOG((1. + zphi_c + zphi_czphi_c)/3.) &
	440	& -1.7320508ATAN((1. + 2.zphi_c)/1.7320508) + 1.813799447
	441	!
	442	zf = zta*zta
	443	zf = zf/(1. + zf)
[11111]	444	zc = MIN(50._wp,0.35_wp*zta)
	445	zstab = 0.5 + SIGN(0.5_wp, zta)
[6723]	446	!
	447	psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)zpsi_k + zfzpsi_c ) &
	448	& - zstab * ( (ABS(1. + 2.zta/3.))*1.5 &
	449	& + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
	450	!
	451	END DO
	452	END DO
	453	!
	454	END FUNCTION psi_h_coare
	455
[11111]	456	!!======================================================================
[6723]	457	END MODULE sbcblk_algo_coare

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

Download in other formats: