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 @ 14017

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

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

Download in other formats: