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_ecmwf.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_ecmwf.F90 @ 14021

Last change on this file since 14021 was 14021, checked in by laurent, 3 years ago
Caught up with trunk rev 14020...
Property svn:keywords set to `Id`
File size: 24.7 KB

Rev	Line
[6723]	1	MODULE sbcblk_algo_ecmwf
	2	!!======================================================================
[12377]	3	!! * MODULE sbcblk_algo_ecmwf *
	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
[13655]	7	!! * the effective bulk wind speed at 10m Ubzu
[6723]	8	!! => all these are used in bulk formulas in sbcblk.F90
	9	!!
[12377]	10	!! Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1)
[6723]	11	!! based on IFS doc (avaible online on the ECMWF's website)
	12	!!
	13	!! Routine turb_ecmwf maintained and developed in AeroBulk
[12377]	14	!! (https://github.com/brodeau/aerobulk)
[6723]	15	!!
[12377]	16	!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
[6723]	17	!!----------------------------------------------------------------------
[6727]	18	!! History : 4.0 ! 2016-02 (L.Brodeau) Original code
[14021]	19	!! 4.2 ! 2020-12 (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
[6723]	20	!!----------------------------------------------------------------------
[6727]	21
	22	!!----------------------------------------------------------------------
[6723]	23	!! turb_ecmwf : computes the bulk turbulent transfer coefficients
	24	!! adjusts t_air and q_air from zt to zu m
	25	!! returns the effective bulk wind speed at 10m
	26	!!----------------------------------------------------------------------
	27	USE dom_oce ! ocean space and time domain
	28	USE phycst ! physical constants
[13655]	29	USE lib_mpp, ONLY: ctl_stop ! distribued memory computing library
	30	USE in_out_manager, ONLY: nit000 ! I/O manager
[13719]	31	USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
[12377]	32	USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB
[6723]	33
	34	IMPLICIT NONE
	35	PRIVATE
	36
[12377]	37	PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF
[6723]	38
[12377]	39	!! ECMWF own values for given constants, taken form IFS documentation...
[13655]	40	REAL(wp), PARAMETER, PUBLIC :: charn0_ecmwf = 0.018_wp ! Charnock constant (pretty high value here !!!
[6723]	41	! ! => Usually 0.011 for moderate winds)
	42	REAL(wp), PARAMETER :: zi0 = 1000. ! scale height of the atmospheric boundary layer...1
	43	REAL(wp), PARAMETER :: Beta0 = 1. ! gustiness parameter ( = 1.25 in COAREv3)
	44	REAL(wp), PARAMETER :: alpha_M = 0.11 ! For roughness length (smooth surface term)
	45	REAL(wp), PARAMETER :: alpha_H = 0.40 ! (Chapter 3, p.34, IFS doc Cy31r1)
	46	REAL(wp), PARAMETER :: alpha_Q = 0.62 !
[12377]	47
[13655]	48	!! * Substitutions
	49	# include "do_loop_substitute.h90"
[12377]	50
[6723]	51	!!----------------------------------------------------------------------
	52	CONTAINS
	53
[12377]	54
	55	SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl)
	56	!!---------------------------------------------------------------------
	57	!! * FUNCTION sbcblk_algo_ecmwf_init *
	58	!!
	59	!! INPUT :
	60	!! -------
	61	!! * l_use_cs : use the cool-skin parameterization
	62	!! * l_use_wl : use the warm-layer parameterization
	63	!!---------------------------------------------------------------------
	64	LOGICAL , INTENT(in) :: l_use_cs ! use the cool-skin parameterization
	65	LOGICAL , INTENT(in) :: l_use_wl ! use the warm-layer parameterization
	66	INTEGER :: ierr
	67	!!---------------------------------------------------------------------
	68	IF( l_use_wl ) THEN
	69	ierr = 0
	70	ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr )
	71	IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' )
	72	dT_wl(:,:) = 0._wp
	73	Hz_wl(:,:) = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars)
	74	ENDIF
	75	IF( l_use_cs ) THEN
	76	ierr = 0
	77	ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
	78	IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' )
	79	dT_cs(:,:) = -0.25_wp ! First guess of skin correction
	80	ENDIF
	81	END SUBROUTINE sbcblk_algo_ecmwf_init
	82
	83
	84
	85	SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
[13655]	86	& Cd, Ch, Ce, t_zu, q_zu, Ubzu, &
	87	& nb_iter, Cdn, Chn, Cen, & ! optional output
[12377]	88	& Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer)
	89	& pdT_wl, pHz_wl ) ! optionals for warm-layer only
[12615]	90	!!----------------------------------------------------------------------------------
[6723]	91	!! * ROUTINE turb_ecmwf *
	92	!!
	93	!! ** Purpose : Computes turbulent transfert coefficients of surface
[12377]	94	!! fluxes according to IFS doc. (cycle 45r1)
[6723]	95	!! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
[12377]	96	!! Returns the effective bulk wind speed at zu to be used in the bulk formulas
[6723]	97	!!
[12377]	98	!! Applies the cool-skin warm-layer correction of the SST to T_s
	99	!! if the net shortwave flux at the surface (Qsw), the downwelling longwave
	100	!! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
	101	!! are provided as (optional) arguments!
[6723]	102	!!
	103	!! INPUT :
	104	!! -------
[12377]	105	!! * kt : current time step (starts at 1)
[6723]	106	!! * zt : height for temperature and spec. hum. of air [m]
[12377]	107	!! * zu : height for wind speed (usually 10m) [m]
[6723]	108	!! * t_zt : potential air temperature at zt [K]
	109	!! * q_zt : specific humidity of air at zt [kg/kg]
[12377]	110	!! * U_zu : scalar wind speed at zu [m/s]
	111	!! * l_use_cs : use the cool-skin parameterization
	112	!! * l_use_wl : use the warm-layer parameterization
[6723]	113	!!
[12377]	114	!! INPUT/OUTPUT:
	115	!! -------------
	116	!! * T_s : always "bulk SST" as input [K]
	117	!! -> unchanged "bulk SST" as output if CSWL not used [K]
	118	!! -> skin temperature as output if CSWL used [K]
[6723]	119	!!
[12377]	120	!! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg]
	121	!! -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True)
	122	!! -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False)
	123	!!
	124	!! OPTIONAL INPUT:
	125	!! ---------------
	126	!! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2]
	127	!! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2]
	128	!! * slp : sea-level pressure [Pa]
	129	!!
	130	!! OPTIONAL OUTPUT:
	131	!! ----------------
	132	!! * pdT_cs : SST increment "dT" for cool-skin correction [K]
	133	!! * pdT_wl : SST increment "dT" for warm-layer correction [K]
	134	!! * pHz_wl : thickness of warm-layer [m]
	135	!!
[6723]	136	!! OUTPUT :
	137	!! --------
	138	!! * Cd : drag coefficient
	139	!! * Ch : sensible heat coefficient
	140	!! * Ce : evaporation coefficient
	141	!! * t_zu : pot. air temperature adjusted at wind height zu [K]
	142	!! * q_zu : specific humidity of air // [kg/kg]
[13655]	143	!! * Ubzu : bulk wind speed at zu [m/s]
[6723]	144	!!
	145	!!
[12377]	146	!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	147	!!----------------------------------------------------------------------------------
[12377]	148	INTEGER, INTENT(in ) :: kt ! current time step
[6723]	149	REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m]
	150	REAL(wp), INTENT(in ) :: zu ! height for U_zu [m]
[12377]	151	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin]
[6723]	152	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin]
[12377]	153	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg]
	154	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg]
[6723]	155	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s]
[12377]	156	LOGICAL , INTENT(in ) :: l_use_cs ! use the cool-skin parameterization
	157	LOGICAL , INTENT(in ) :: l_use_wl ! use the warm-layer parameterization
[6723]	158	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau)
	159	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens)
	160	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat)
	161	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K]
	162	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg]
[13655]	163	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s]
[6723]	164	!
[13655]	165	INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations
	166	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN
	167	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN
	168	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN
[12377]	169	REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2]
	170	REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2]
	171	REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa]
	172	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs
	173	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K]
	174	REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m]
	175	!
[13655]	176	INTEGER :: nbit, jit
[12377]	177	LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U
[6723]	178	!
[12615]	179	REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star
	180	REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu
	181	REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air
[12377]	182	REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length...
	183	REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q
[6723]	184	!
[12377]	185	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST
[6723]	186	!
[12377]	187	REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h
	188	REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
	189	CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
	190	!!----------------------------------------------------------------------------------
	191	IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl)
	192
[13655]	193	nbit = nb_iter0
	194	IF( PRESENT(nb_iter) ) nbit = nb_iter
	195
[12615]	196	l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
[6723]	197
[12377]	198	!! Initializations for cool skin and warm layer:
	199	IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
	200	& CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
[6723]	201
[12377]	202	IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
	203	& CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
	204
	205	IF( l_use_cs .OR. l_use_wl ) THEN
	206	ALLOCATE ( zsst(jpi,jpj) )
	207	zsst = T_s ! backing up the bulk SST
	208	IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction
	209	q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
	210	ENDIF
	211
	212
	213	! Identical first gess as in COARE, with IFS parameter values though...
	214	!
[6723]	215	!! First guess of temperature and humidity at height zu:
[12377]	216	t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions...
	217	q_zu = MAX( q_zt , 1.e-6_wp ) ! "
[6723]	218
	219	!! Pot. temp. difference (and we don't want it to be 0!)
[12377]	220	dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
	221	dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
[6723]	222
[12377]	223	znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
[6723]	224
[13655]	225	Ubzu = SQRT(U_zuU_zu + 0.5_wp0.5_wp) ! initial guess for wind gustiness contribution
[6723]	226
[12377]	227	ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
	228	ztmp1 = LOG(10._wp*10000._wp) ! " " "
[13655]	229	u_star = 0.035_wpUbzuztmp1/ztmp0 ! (u* = 0.035*Un10)
[6723]	230
[13655]	231	z0 = charn0_ecmwfu_staru_star/grav + 0.11_wp*znu_a/u_star
[12377]	232	z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on)
[6723]	233
[12377]	234	z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
	235	z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on)
	236
[13655]	237	Cd = MAX( (vkarmn/ztmp0)**2 , Cx_min ) ! first guess of Cd
[6723]	238
[13655]	239	ztmp0 = vkarmn2/LOG(zt/z0t)/Cd
[6723]	240
[13655]	241	ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN)
[6723]	242
[12377]	243	!! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
	244	ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
[13655]	245	func_h = (1._wp - ztmp1) * ztmp0ztmp2 / (1._wp - ztmp2zi00.004_wpBeta0**3/zu) & ! BRN < 0
	246	& + ztmp1 * ( ztmp0ztmp2 + 27._wp/9._wpztmp2*ztmp2 ) ! BRN > 0
[14017]	247
[6723]	248	!! First guess M-O stability dependent scaling params.(u,t,q*) to estimate z0 and z/L
[12377]	249	ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
[6723]	250
[13655]	251	u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on)
[6723]	252	t_star = dt_zu*ztmp0
	253	q_star = dq_zu*ztmp0
	254
[12377]	255	! What needs to be done if zt /= zu:
[6723]	256	IF( .NOT. l_zt_equal_zu ) THEN
	257	!! First update of values at zu (or zt for wind)
	258	ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(ztfunc_h/zu) ! ztfunc_h/zu == zeta_t
[12377]	259	ztmp1 = LOG(zt/zu) + ztmp0
[6723]	260	t_zu = t_zt - t_star/vkarmn*ztmp1
	261	q_zu = q_zt - q_star/vkarmn*ztmp1
[12377]	262	q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
[6723]	263	!
[12377]	264	dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
	265	dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
[6723]	266	ENDIF
	267
	268
	269	!! => that was same first guess as in COARE...
	270
	271
[13655]	272	!! First guess of inverse of Obukov length (1/L) :
[12377]	273	Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
[6723]	274
[13655]	275	!! Functions such as u* = Ubzu*vkarmn/func_m
[12377]	276	ztmp0 = zu*Linv
	277	func_m = LOG(zu) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv)
	278	func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
[6723]	279
	280	!! ITERATION BLOCK
[13655]	281	DO jit = 1, nbit
[6723]	282
	283	!! Bulk Richardson Number at z=zu (Eq. 3.25)
[13655]	284	ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN)
[6723]	285
[13655]	286	!! New estimate of the inverse of the Obukhon length (Linv == zeta/zu) :
[12377]	287	Linv = ztmp0func_mfunc_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
	288	!! Note: it is slightly different that the L we would get with the usual
	289	Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...)
[6723]	290
	291	!! Update func_m with new Linv:
[12377]	292	func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zuLinv) + psi_m_ecmwf(z0Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu!
[6723]	293
	294	!! Need to update roughness lengthes:
[13655]	295	u_star = Ubzu*vkarmn/func_m
[6723]	296	ztmp2 = u_star*u_star
	297	ztmp1 = znu_a/u_star
[13655]	298	z0 = MIN( ABS( alpha_Mztmp1 + charn0_ecmwfztmp2/grav ) , 0.001_wp)
	299	z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
	300	z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp)
[6723]	301
[12377]	302	!! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8)
	303	ztmp2 = Beta0Beta0ztmp2(MAX(-zi0Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
	304	!! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
[13655]	305	Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed
	306	! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0.
[6723]	307
	308
	309	!! Need to update "theta" and "q" at zu in case they are given at different heights
	310	!! as well the air-sea differences:
	311	IF( .NOT. l_zt_equal_zu ) THEN
	312	!! Arrays func_m and func_h are free for a while so using them as temporary arrays...
[12377]	313	func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!!
	314	func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!!
[6723]	315
	316	ztmp2 = psi_h_ecmwf(z0t*Linv)
	317	ztmp0 = func_h - ztmp2
[12377]	318	ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
[6723]	319	t_star = dt_zu*ztmp1
	320	ztmp2 = ztmp0 - func_m + ztmp2
	321	ztmp1 = LOG(zt/zu) + ztmp2
	322	t_zu = t_zt - t_star/vkarmn*ztmp1
	323
	324	ztmp2 = psi_h_ecmwf(z0q*Linv)
	325	ztmp0 = func_h - ztmp2
[12377]	326	ztmp1 = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0)
[6723]	327	q_star = dq_zu*ztmp1
	328	ztmp2 = ztmp0 - func_m + ztmp2
[12377]	329	ztmp1 = LOG(zt/zu) + ztmp2
[6723]	330	q_zu = q_zt - q_star/vkarmn*ztmp1
[12377]	331	ENDIF
[6723]	332
[12377]	333	!! Updating because of updated z0 and z0t and new Linv...
	334	ztmp0 = zu*Linv
	335	func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
	336	func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
[9019]	337
[6723]	338
[12377]	339	IF( l_use_cs ) THEN
	340	!! Cool-skin contribution
[6723]	341
[13655]	342	CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
[12377]	343	& ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0
[6723]	344
[12377]	345	CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst ) ! Qnsol -> ztmp1
	346
	347	T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1)
	348	IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1)
	349	q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
	350
	351	ENDIF
	352
	353	IF( l_use_wl ) THEN
	354	!! Warm-layer contribution
[13655]	355	CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
[12377]	356	& ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2
	357	CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst )
	358	!! Updating T_s and q_s !!!
	359	T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) !
	360	IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
	361	q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
	362	ENDIF
	363
	364	IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
	365	dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
	366	dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
	367	ENDIF
	368
[13655]	369	END DO !DO jit = 1, nbit
[12377]	370
[13655]	371	Cd = MAX( vkarmn2/(func_m*func_m) , Cx_min )
	372	Ch = MAX( vkarmn2/(func_m*func_h) , Cx_min )
	373	ztmp2 = LOG(zu/z0q) - psi_h_ecmwf(zuLinv) + psi_h_ecmwf(z0qLinv) ! func_q
	374	Ce = MAX( vkarmn2/(func_m*ztmp2) , Cx_min )
[6723]	375
[13655]	376	IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min )
	377	IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min )
	378	IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0q)*LOG(zu/z0q)) , Cx_min )
[9019]	379
[12377]	380	IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
	381	IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
	382	IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
[6723]	383
[12377]	384	IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
[6723]	385
[12377]	386	END SUBROUTINE turb_ecmwf
	387
	388
[6723]	389	FUNCTION psi_m_ecmwf( pzeta )
	390	!!----------------------------------------------------------------------------------
	391	!! Universal profile stability function for momentum
	392	!! ECMWF / as in IFS cy31r1 documentation, available online
	393	!! at ecmwf.int
	394	!!
	395	!! pzeta : stability paramenter, z/L where z is altitude measurement
	396	!! and L is M-O length
	397	!!
[12377]	398	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	399	!!----------------------------------------------------------------------------------
	400	REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf
	401	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
	402	!
	403	INTEGER :: ji, jj ! dummy loop indices
[13655]	404	REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc
[6723]	405	!!----------------------------------------------------------------------------------
[13655]	406	zc = 5._wp/0.35_wp
[13460]	407	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13655]	408	!
	409	zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
	410
	411	! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
	412	zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5
	413	zx = SQRT(zx2) ! (1 - 16z)^0.25
	414	ztmp = 1._wp + zx
	415	zpsi_unst = LOG( 0.125_wpztmpztmp(1._wp + zx2) ) - 2._wpATAN( zx ) + 0.5_wp*rpi
	416
	417	! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
	418	zpsi_stab = -2._wp/3._wp(zta - zc)EXP(-0.35_wp*zta) &
	419	& - zta - 2._wp/3._wp*zc
	420	!
	421	zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
	422	!
	423	psi_m_ecmwf(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
	424	& + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
	425	!
[12377]	426	END_2D
[6723]	427	END FUNCTION psi_m_ecmwf
	428
[12377]	429
[6723]	430	FUNCTION psi_h_ecmwf( pzeta )
	431	!!----------------------------------------------------------------------------------
	432	!! Universal profile stability function for temperature and humidity
	433	!! ECMWF / as in IFS cy31r1 documentation, available online
	434	!! at ecmwf.int
	435	!!
	436	!! pzeta : stability paramenter, z/L where z is altitude measurement
	437	!! and L is M-O length
	438	!!
[12377]	439	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]	440	!!----------------------------------------------------------------------------------
	441	REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
	442	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
	443	!
	444	INTEGER :: ji, jj ! dummy loop indices
[13655]	445	REAL(wp) :: zta, zx2, zpsi_unst, zpsi_stab, zstab, zc
[6723]	446	!!----------------------------------------------------------------------------------
[13655]	447	zc = 5._wp/0.35_wp
[6723]	448	!
[13460]	449	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[13655]	450	!
	451	zta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!):
	452	!
	453	! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
	454	zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5
	455	zpsi_unst = 2._wpLOG( 0.5_wp(1._wp + zx2) )
	456	!
	457	! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
	458	zpsi_stab = -2._wp/3._wp(zta - zc)EXP(-0.35_wp*zta) &
	459	& - ABS(1._wp + 2._wp/3._wpzta)1.5_wp - 2._wp/3._wpzc + 1._wp
	460	!
	461	! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
	462	!
	463	zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
	464	!
	465	psi_h_ecmwf(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable
[14017]	466	& + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable
[13655]	467	!
[12377]	468	END_2D
[6723]	469	END FUNCTION psi_h_ecmwf
	470
	471
	472	!!======================================================================
	473	END MODULE sbcblk_algo_ecmwf

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

Download in other formats: