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, 14 months ago

LB: gooodbye "coare3p5", hello "coare3p6" + cleaner air humidity handling + support for Relative humidity.

File size: 18.5 KB
Line 
1MODULE sbcblk_phy
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_phy  ***
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
10   !!   virt_temp     : virtual (aka sensible) temperature (potential or absolute)
11   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
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
14   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
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
18   !!   q_sat         : saturation humidity as a function of SLP and temperature
19   !!   q_air_rh      : specific humidity as a function of RH, t_air and SLP
20   
21   USE dom_oce        ! ocean space and time domain
22   USE phycst         ! physical constants
23
24   IMPLICIT NONE
25   PRIVATE
26
27   INTERFACE gamma_moist
28      MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr
29   END INTERFACE gamma_moist
30
31   PUBLIC virt_temp
32   PUBLIC rho_air
33   PUBLIC visc_air
34   PUBLIC L_vap
35   PUBLIC cp_air
36   PUBLIC gamma_moist
37   PUBLIC One_on_L
38   PUBLIC Ri_bulk
39   PUBLIC q_sat
40   PUBLIC q_air_rh
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   !!----------------------------------------------------------------------
47CONTAINS
48
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
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      !!
80      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
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_dry*ptak * ( 1._wp + rctv0*pqa )  )
89      !
90   END FUNCTION rho_air
91
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-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
110         END DO
111      END DO
112      !
113   END FUNCTION visc_air
114
115   FUNCTION L_vap( psst )
116      !!---------------------------------------------------------------------------------
117      !!                           ***  FUNCTION L_vap  ***
118      !!
119      !! ** Purpose : Compute the latent heat of vaporization of water out of temperature
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
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
147   FUNCTION gamma_moist_vctr( ptak, pqa )
148      !!----------------------------------------------------------------------------------
149      !!                           ***  FUNCTION gamma_moist_vctr  ***
150      !!
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
154      !!
155      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
156      !!----------------------------------------------------------------------------------
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
160      !
161      INTEGER  ::   ji, jj         ! dummy loop indices
162      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
163      !!----------------------------------------------------------------------------------
164      !
165      DO jj = 1, jpj
166         DO ji = 1, jpi
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) !    "                   "                     "
169            !
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 + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
173         END DO
174      END DO
175      !
176   END FUNCTION gamma_moist_vctr
177
178   FUNCTION gamma_moist_sclr( ptak, pqa )
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      !!
184      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
185      !!----------------------------------------------------------------------------------
186      REAL(wp)             :: gamma_moist_sclr
187      REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
188      !
189      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
190      !!----------------------------------------------------------------------------------
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 + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/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
213      !
214      INTEGER  ::   ji, jj         ! dummy loop indices
215      REAL(wp) ::     zqa          ! local scalar
216      !!-------------------------------------------------------------------
217      !
218      DO jj = 1, jpj
219         DO ji = 1, jpi
220            !
221            zqa = (1._wp + rctv0*pqa(ji,jj))
222            !
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) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
229               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
230            !
231         END DO
232      END DO
233      !
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
237
238   FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub )
239      !!----------------------------------------------------------------------------------
240      !! Bulk Richardson number according to "wide-spread equation"...
241      !!
242      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
243      !!----------------------------------------------------------------------------------
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]
251      !
252      INTEGER  ::   ji, jj                                ! dummy loop indices
253      REAL(wp) ::   zqa, zta, zgamma, zdth_v, ztv, zsstv  ! local scalars
254      !!-------------------------------------------------------------------
255      !
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 + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!)
265            !
266            zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature"
267            !
268            ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer
269            !
270            Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk
271            !
272         END DO
273      END DO
274   END FUNCTION Ri_bulk
275
276
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.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  &
301         &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
302      !
303   END FUNCTION e_sat_sclr
304
305
306   FUNCTION q_sat( ptak, pslp )
307      !!----------------------------------------------------------------------------------
308      !!                           ***  FUNCTION q_sat  ***
309      !!
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      !!
314      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
315      !!----------------------------------------------------------------------------------
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]
319      !
320      INTEGER  ::   ji, jj         ! dummy loop indices
321      REAL(wp) ::   ze_sat   ! local scalar
322      !!----------------------------------------------------------------------------------
323      !
324      DO jj = 1, jpj
325         DO ji = 1, jpi
326            !
327            ze_sat =  e_sat_sclr( ptak(ji,jj), pslp(ji,jj) )
328            !
329            q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat )
330            !
331         END DO
332      END DO
333      !
334   END FUNCTION q_sat
335
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) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze)
355         END DO
356      END DO
357      !
358   END FUNCTION q_air_rh
359
360
361   !!======================================================================
362END MODULE sbcblk_phy
Note: See TracBrowser for help on using the repository browser.