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

Context Navigation

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

Last change on this file since 11284 was 11284, checked in by laurent, 5 years ago
LB: gooodbye "coare3p5", hello "coare3p6" + cleaner air humidity handling + support for Relative humidity.
File size: 18.5 KB

Rev	Line
[11182]	1	MODULE sbcblk_phy
[11111]	2	!!======================================================================
[11182]	3	!! * MODULE sbcblk_phy *
[11111]	4	!! A set of functions to compute air themodynamics parameters
	5	!! needed by Aerodynamic Bulk Formulas
	6	!!=====================================================================
	7	!! 4.0 ! 2019 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
	8	!!----------------------------------------------------------------------
	9
[11209]	10	!! virt_temp : virtual (aka sensible) temperature (potential or absolute)
[11111]	11	!! rho_air : density of (moist) air (depends on T_air, q_air and SLP
[11209]	12	!! visc_air : kinematic viscosity (aka Nu_air) of air from temperature
	13	!! L_vap : latent heat of vaporization of water as a function of temperature
[11111]	14	!! cp_air : specific heat of (moist) air (depends spec. hum. q_air)
[11209]	15	!! gamma_moist : adiabatic lapse-rate of moist air
	16	!! One_on_L : 1. / ( Monin-Obukhov length )
	17	!! Ri_bulk : bulk Richardson number aka BRN
[11111]	18	!! q_sat : saturation humidity as a function of SLP and temperature
[11284]	19	!! q_air_rh : specific humidity as a function of RH, t_air and SLP
	20
[11111]	21	USE dom_oce ! ocean space and time domain
	22	USE phycst ! physical constants
	23
	24	IMPLICIT NONE
	25	PRIVATE
	26
[11209]	27	INTERFACE gamma_moist
	28	MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr
	29	END INTERFACE gamma_moist
	30
	31	PUBLIC virt_temp
[11111]	32	PUBLIC rho_air
[11209]	33	PUBLIC visc_air
	34	PUBLIC L_vap
[11111]	35	PUBLIC cp_air
[11209]	36	PUBLIC gamma_moist
	37	PUBLIC One_on_L
	38	PUBLIC Ri_bulk
[11111]	39	PUBLIC q_sat
[11284]	40	PUBLIC q_air_rh
[11111]	41
	42	!!----------------------------------------------------------------------
	43	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
	44	!! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
	45	!! Software governed by the CeCILL license (see ./LICENSE)
	46	!!----------------------------------------------------------------------
	47	CONTAINS
	48
[11209]	49	FUNCTION virt_temp( pta, pqa )
	50	!!------------------------------------------------------------------------
	51	!!
	52	!! Compute the (absolute/potential) virtual temperature, knowing the
	53	!! (absolute/potential) temperature and specific humidity
	54	!!
	55	!! If input temperature is absolute then output vitual temperature is absolute
	56	!! If input temperature is potential then output vitual temperature is potential
	57	!!
	58	!! Author: L. Brodeau, June 2019 / AeroBulk
	59	!! (https://github.com/brodeau/aerobulk/)
	60	!!------------------------------------------------------------------------
	61	REAL(wp), DIMENSION(jpi,jpj) :: virt_temp !: 1./(Monin Obukhov length) [m^-1]
	62	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta, & !: absolute or potetntial air temperature [K]
	63	& pqa !: specific humidity of air [kg/kg]
	64	!!-------------------------------------------------------------------
	65	!
	66	virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
	67	!!
	68	!! This is exactly the same sing that:
	69	!! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa))
	70	!! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
	71	!
	72	END FUNCTION virt_temp
	73
[11111]	74	FUNCTION rho_air( ptak, pqa, pslp )
	75	!!-------------------------------------------------------------------------------
	76	!! * FUNCTION rho_air *
	77	!!
	78	!! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
	79	!!
[11209]	80	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[11111]	81	!!-------------------------------------------------------------------------------
	82	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
	83	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
	84	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]
	85	REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]
	86	!!-------------------------------------------------------------------------------
	87	!
	88	rho_air = pslp / ( R_dryptak ( 1._wp + rctv0*pqa ) )
	89	!
	90	END FUNCTION rho_air
	91
[11209]	92	FUNCTION visc_air(ptak)
	93	!!----------------------------------------------------------------------------------
	94	!! Air kinetic viscosity (m^2/s) given from temperature in degrees...
	95	!!
	96	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
	97	!!----------------------------------------------------------------------------------
	98	REAL(wp), DIMENSION(jpi,jpj) :: visc_air !
	99	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)
	100	!
	101	INTEGER :: ji, jj ! dummy loop indices
	102	REAL(wp) :: ztc, ztc2 ! local scalar
	103	!!----------------------------------------------------------------------------------
	104	!
	105	DO jj = 1, jpj
	106	DO ji = 1, jpi
	107	ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C
	108	ztc2 = ztc*ztc
	109	visc_air(ji,jj) = 1.326e-5(1. + 6.542E-3ztc + 8.301e-6ztc2 - 4.84e-9ztc2*ztc)
	110	END DO
	111	END DO
	112	!
	113	END FUNCTION visc_air
[11111]	114
[11209]	115	FUNCTION L_vap( psst )
	116	!!---------------------------------------------------------------------------------
	117	!! * FUNCTION L_vap *
	118	!!
[11284]	119	!! ** Purpose : Compute the latent heat of vaporization of water out of temperature
[11209]	120	!!
	121	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
	122	!!----------------------------------------------------------------------------------
	123	REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]
	124	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
	125	!!----------------------------------------------------------------------------------
	126	!
	127	L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6
	128	!
	129	END FUNCTION L_vap
	130
[11111]	131	FUNCTION cp_air( pqa )
	132	!!-------------------------------------------------------------------------------
	133	!! * FUNCTION cp_air *
	134	!!
	135	!! ** Purpose : Compute specific heat (Cp) of moist air
	136	!!
	137	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
	138	!!-------------------------------------------------------------------------------
	139	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
	140	REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg]
	141	!!-------------------------------------------------------------------------------
	142	!
	143	cp_air = rCp_dry + rCp_vap * pqa
	144	!
	145	END FUNCTION cp_air
	146
[11209]	147	FUNCTION gamma_moist_vctr( ptak, pqa )
[11111]	148	!!----------------------------------------------------------------------------------
[11209]	149	!! * FUNCTION gamma_moist_vctr *
[11111]	150	!!
[11209]	151	!! ** Purpose : Compute the moist adiabatic lapse-rate.
	152	!! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
	153	!! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
[11111]	154	!!
	155	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
	156	!!----------------------------------------------------------------------------------
[11209]	157	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
	158	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]
	159	REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist_vctr ! moist adiabatic lapse-rate
[11111]	160	!
	161	INTEGER :: ji, jj ! dummy loop indices
[11209]	162	REAL(wp) :: zta, zqa, zwa, ziRT ! local scalar
[11111]	163	!!----------------------------------------------------------------------------------
	164	!
	165	DO jj = 1, jpj
	166	DO ji = 1, jpi
[11209]	167	zta = MAX( ptak(ji,jj), 180._wp) ! prevents screw-up over masked regions where field == 0.
	168	zqa = MAX( pqa(ji,jj), 1.E-6_wp) ! " " "
[11111]	169	!
[11209]	170	zwa = zqa / (1. - zqa) ! w is mixing ratio w = q/(1-q) \| q = w/(1+w)
	171	ziRT = 1._wp/(R_dry*zta) ! 1/RT
	172	gamma_moist_vctr(ji,jj) = grav * ( 1._wp + rLevapzwaziRT ) / ( rCp_dry + rLevaprLevapzwareps0ziRT/zta )
[11111]	173	END DO
	174	END DO
	175	!
[11209]	176	END FUNCTION gamma_moist_vctr
[11111]	177
[11209]	178	FUNCTION gamma_moist_sclr( ptak, pqa )
[11111]	179	!!----------------------------------------------------------------------------------
	180	!! ** Purpose : Compute the moist adiabatic lapse-rate.
	181	!! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
	182	!! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
	183	!!
[11209]	184	!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[11111]	185	!!----------------------------------------------------------------------------------
[11209]	186	REAL(wp) :: gamma_moist_sclr
	187	REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
[11111]	188	!
[11209]	189	REAL(wp) :: zta, zqa, zwa, ziRT ! local scalar
[11111]	190	!!----------------------------------------------------------------------------------
[11209]	191	zta = MAX( ptak, 180._wp) ! prevents screw-up over masked regions where field == 0.
	192	zqa = MAX( pqa, 1.E-6_wp) ! " " "
	193	!!
	194	zwa = zqa / (1._wp - zqa) ! w is mixing ratio w = q/(1-q) \| q = w/(1+w)
	195	ziRT = 1._wp / (R_dry*zta) ! 1/RT
	196	gamma_moist_sclr = grav * ( 1._wp + rLevapzwaziRT ) / ( rCp_dry + rLevaprLevapzwareps0ziRT/zta )
	197	!!
	198	END FUNCTION gamma_moist_sclr
	199
	200	FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
	201	!!------------------------------------------------------------------------
	202	!!
	203	!! Evaluates the 1./(Monin Obukhov length) from air temperature and
	204	!! specific humidity, and frictional scales u, t and q*
	205	!!
	206	!! Author: L. Brodeau, June 2016 / AeroBulk
	207	!! (https://github.com/brodeau/aerobulk/)
	208	!!------------------------------------------------------------------------
	209	REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1]
	210	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha, & !: average potetntial air temperature [K]
	211	& pqa, & !: average specific humidity of air [kg/kg]
	212	& pus, pts, pqs !: frictional velocity, temperature and humidity
[11111]	213	!
[11209]	214	INTEGER :: ji, jj ! dummy loop indices
	215	REAL(wp) :: zqa ! local scalar
	216	!!-------------------------------------------------------------------
	217	!
[11111]	218	DO jj = 1, jpj
	219	DO ji = 1, jpi
[11209]	220	!
	221	zqa = (1._wp + rctv0*pqa(ji,jj))
	222	!
[11266]	223	! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
	224	! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
	225	! or
	226	! b/ -u* [ theta* + 0.61 theta q* ]
	227	!
	228	One_on_L(ji,jj) = gravvkarmn( pts(ji,jj)zqa + rctv0ptha(ji,jj)*pqs(ji,jj) ) &
[11209]	229	& / MAX( pus(ji,jj)pus(ji,jj)ptha(ji,jj)*zqa , 1.E-9_wp )
	230	!
[11111]	231	END DO
	232	END DO
	233	!
[11209]	234	One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
	235	!
	236	END FUNCTION One_on_L
[11111]	237
[11209]	238	FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub )
	239	!!----------------------------------------------------------------------------------
	240	!! Bulk Richardson number according to "wide-spread equation"...
[11111]	241	!!
[11209]	242	!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
[11111]	243	!!----------------------------------------------------------------------------------
[11209]	244	REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk
	245	REAL(wp) , INTENT(in) :: pz ! height above the sea (aka "delta z") [m]
	246	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! SST [K]
	247	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K]
	248	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg]
	249	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air spec. hum. at height "pz" [kg/kg]
	250	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub ! bulk wind speed [m/s]
[11111]	251	!
[11209]	252	INTEGER :: ji, jj ! dummy loop indices
	253	REAL(wp) :: zqa, zta, zgamma, zdth_v, ztv, zsstv ! local scalars
	254	!!-------------------------------------------------------------------
[11111]	255	!
[11209]	256	DO jj = 1, jpj
	257	DO ji = 1, jpi
	258	!
	259	zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj)) ! ~ mean q within the layer...
	260	zta = 0.5_wp( psst(ji,jj) + ptha(ji,jj) - gamma_moist(ptha(ji,jj),zqa)pz ) ! ~ mean absolute temperature of air within the layer
	261	zta = 0.5_wp( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta, zqa)pz ) ! ~ mean absolute temperature of air within the layer
	262	zgamma = gamma_moist(zta, zqa) ! Adiabatic lapse-rate for moist air within the layer
	263	!
	264	zsstv = psst(ji,jj)(1._wp + rctv0pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!)
	265	!
	266	zdth_v = ptha(ji,jj)(1._wp + rctv0pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature"
	267	!
	268	ztv = 0.5_wp( zsstv + (ptha(ji,jj) - zgammapz)(1._wp + rctv0pqa(ji,jj)) ) ! ~ mean absolute virtual temp. within the layer
	269	!
	270	Ri_bulk(ji,jj) = gravzdth_vpz / ( ztvpub(ji,jj)pub(ji,jj) ) ! the usual definition of Ri_bulk
	271	!
	272	END DO
	273	END DO
	274	END FUNCTION Ri_bulk
[11111]	275
	276
[11284]	277	FUNCTION e_sat_sclr( ptak, pslp )
	278	!!----------------------------------------------------------------------------------
	279	!! * FUNCTION e_sat_sclr *
	280	!! < SCALAR argument version >
	281	!! ** Purpose : water vapor at saturation in [Pa]
	282	!! Based on accurate estimate by Goff, 1957
	283	!!
	284	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
	285	!!
	286	!! Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here
	287	!!----------------------------------------------------------------------------------
	288	REAL(wp), INTENT(in) :: ptak ! air temperature [K]
	289	REAL(wp), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]
	290	REAL(wp) :: e_sat_sclr ! water vapor at saturation [kg/kg]
	291	!
	292	REAL(wp) :: zta, ztmp ! local scalar
	293	!!----------------------------------------------------------------------------------
	294	!
	295	zta = MAX( ptak , 180._wp ) ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
	296	ztmp = rt0 / zta
	297	!
	298	! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957)
	299	e_sat_sclr = 100.( 10.( 10.79574(1. - ztmp) - 5.028*LOG10(zta/rt0) &
	300	& + 1.5047510.(-4)(1. - 10.*(-8.2969(zta/rt0 - 1.)) ) &
	301	& + 0.4287310.(-3)(10.*(4.76955(1. - ztmp)) - 1.) + 0.78614) )
	302	!
	303	END FUNCTION e_sat_sclr
	304
	305
[11209]	306	FUNCTION q_sat( ptak, pslp )
[11111]	307	!!----------------------------------------------------------------------------------
[11209]	308	!! * FUNCTION q_sat *
[11111]	309	!!
[11209]	310	!! ** Purpose : Specific humidity at saturation in [kg/kg]
	311	!! Based on accurate estimate of "e_sat"
	312	!! aka saturation water vapor (Goff, 1957)
	313	!!
[11111]	314	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
	315	!!----------------------------------------------------------------------------------
[11209]	316	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
	317	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]
	318	REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]
[11111]	319	!
[11209]	320	INTEGER :: ji, jj ! dummy loop indices
[11284]	321	REAL(wp) :: ze_sat ! local scalar
[11111]	322	!!----------------------------------------------------------------------------------
	323	!
	324	DO jj = 1, jpj
	325	DO ji = 1, jpi
[11209]	326	!
[11284]	327	ze_sat = e_sat_sclr( ptak(ji,jj), pslp(ji,jj) )
[11209]	328	!
[11284]	329	q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat )
[11209]	330	!
[11111]	331	END DO
	332	END DO
	333	!
[11209]	334	END FUNCTION q_sat
[11111]	335
[11284]	336	FUNCTION q_air_rh(prha, ptak, pslp)
	337	!!----------------------------------------------------------------------------------
	338	!! Specific humidity of air out of Relative Humidity
	339	!!
	340	!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
	341	!!----------------------------------------------------------------------------------
	342	REAL(wp), DIMENSION(jpi,jpj) :: q_air_rh
	343	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha !: relative humidity [fraction, not %!!!]
	344	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak !: air temperature [K]
	345	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp !: atmospheric pressure [Pa]
	346	!
	347	INTEGER :: ji, jj ! dummy loop indices
	348	REAL(wp) :: ze ! local scalar
	349	!!----------------------------------------------------------------------------------
	350	!
	351	DO jj = 1, jpj
	352	DO ji = 1, jpi
	353	ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj))
	354	q_air_rh(ji,jj) = zereps0/(pslp(ji,jj) - (1. - reps0)ze)
	355	END DO
	356	END DO
	357	!
	358	END FUNCTION q_air_rh
[11209]	359
[11284]	360
[11111]	361	!!======================================================================
[11182]	362	END MODULE sbcblk_phy

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

Download in other formats: