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

Context Navigation

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_algo_ecmwf.F90 @ 12154

Last change on this file since 12154 was 12154, checked in by cetlod, 4 years ago
commit
Property svn:keywords set to `Id`
File size: 24.9 KB

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

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

Download in other formats: