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/tickets_icb_1900/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcblk_algo_ecmwf.F90 @ 14016

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

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

Download in other formats: