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

Last change on this file since 11845 was 11845, checked in by laurent, 11 months ago

Improving syntax consistency

File size: 35.3 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 (fraction, not %), 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   !!  (mainly removed from sbcblk.F90)
28   REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp   !: Specic heat of dry air, constant pressure      [J/K/kg]
29   REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp   !: Specic heat of water vapor, constant pressure  [J/K/kg]
30   REAL(wp), PARAMETER, PUBLIC :: R_dry   = 287.05_wp   !: Specific gas constant for dry air              [J/K/kg]
31   REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg]
32   REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
33   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
34   REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...)
35   REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp   !: transfer coefficient over ice
36   REAL(wp), PARAMETER, PUBLIC :: albo    = 0.066_wp    !: ocean albedo assumed to be constant
37   !
38   REAL(wp), PARAMETER, PUBLIC :: rho0_a  = 1.2_wp      !: Approx. of density of air                       [kg/m^3]
39   REAL(wp), PARAMETER, PUBLIC :: rho0_w  = 1025._wp    !: Density of sea-water  (ECMWF->1025)             [kg/m^3]
40   REAL(wp), PARAMETER, PUBLIC :: radrw = rho0_a/rho0_w !: Density ratio
41   REAL(wp), PARAMETER, PUBLIC :: sq_radrw = SQRT(rho0_a/rho0_w)
42   REAL(wp), PARAMETER, PUBLIC :: rCp0_w  = 4190._wp    !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg]
43   REAL(wp), PARAMETER, PUBLIC :: rnu0_w  = 1.e-6_wp    !: kinetic viscosity of water                      [m^2/s]
44   REAL(wp), PARAMETER, PUBLIC :: rk0_w   = 0.6_wp      !: thermal conductivity of water (at 20C)          [W/m/K]
45   !
46   REAL(wp), PARAMETER, PUBLIC :: emiss_w = 1._wp       !: Surface emissivity (black-body long-wave radiation) of sea-water []
47   !                                                    !: Theoretically close to 0.97! Yet, taken equal as 1 to account for
48   !                                                    !: the small fraction of downwelling shortwave reflected at the
49   !                                                    !: surface (Lind & Katsaros, 1986)
50   REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp  !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt
51   REAL(wp), PARAMETER, PUBLIC :: rtt0 = 273.16_wp        !: triple point of temperature    [K]
52   !
53   REAL(wp), PARAMETER, PUBLIC :: rcst_cs = -16._wp*9.80665_wp*rho0_w*rCp0_w*rnu0_w*rnu0_w*rnu0_w/(rk0_w*rk0_w) !: for cool-skin parameterizations... (grav = 9.80665_wp)
54   !                              => see eq.(14) in Fairall et al. 1996   (eq.(6) of Zeng aand Beljaars is WRONG! (typo?)
55
56
57   INTERFACE gamma_moist
58      MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr
59   END INTERFACE gamma_moist
60
61   INTERFACE e_sat
62      MODULE PROCEDURE e_sat_vctr, e_sat_sclr
63   END INTERFACE e_sat
64
65   INTERFACE L_vap
66      MODULE PROCEDURE L_vap_vctr, L_vap_sclr
67   END INTERFACE L_vap
68
69   INTERFACE rho_air
70      MODULE PROCEDURE rho_air_vctr, rho_air_sclr
71   END INTERFACE rho_air
72
73   INTERFACE cp_air
74      MODULE PROCEDURE cp_air_vctr, cp_air_sclr
75   END INTERFACE cp_air
76
77   INTERFACE alpha_sw
78      MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr
79   END INTERFACE alpha_sw
80
81   INTERFACE bulk_formula
82      MODULE PROCEDURE bulk_formula_vctr, bulk_formula_sclr
83   END INTERFACE bulk_formula
84
85
86
87   PUBLIC virt_temp
88   PUBLIC rho_air
89   PUBLIC visc_air
90   PUBLIC L_vap
91   PUBLIC cp_air
92   PUBLIC gamma_moist
93   PUBLIC One_on_L
94   PUBLIC Ri_bulk
95   PUBLIC q_sat
96   PUBLIC q_air_rh
97   !:
98   PUBLIC update_qnsol_tau
99   PUBLIC alpha_sw
100   PUBLIC bulk_formula
101
102   !!----------------------------------------------------------------------
103   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
104   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
105   !! Software governed by the CeCILL license (see ./LICENSE)
106   !!----------------------------------------------------------------------
107CONTAINS
108
109   FUNCTION virt_temp( pta, pqa )
110      !!------------------------------------------------------------------------
111      !!
112      !! Compute the (absolute/potential) virtual temperature, knowing the
113      !! (absolute/potential) temperature and specific humidity
114      !!
115      !! If input temperature is absolute then output vitual temperature is absolute
116      !! If input temperature is potential then output vitual temperature is potential
117      !!
118      !! Author: L. Brodeau, June 2019 / AeroBulk
119      !!         (https://github.com/brodeau/aerobulk/)
120      !!------------------------------------------------------------------------
121      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp         !: 1./(Monin Obukhov length) [m^-1]
122      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta,  &  !: absolute or potetntial air temperature [K]
123         &                                        pqa      !: specific humidity of air   [kg/kg]
124      !!-------------------------------------------------------------------
125      !
126      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
127      !!
128      !! This is exactly the same sing that:
129      !! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa))
130      !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
131      !
132   END FUNCTION virt_temp
133
134   FUNCTION rho_air_vctr( ptak, pqa, pslp )
135      !!-------------------------------------------------------------------------------
136      !!                           ***  FUNCTION rho_air_vctr  ***
137      !!
138      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
139      !!
140      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
141      !!-------------------------------------------------------------------------------
142      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
143      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
144      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
145      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air_vctr   ! density of moist air   [kg/m^3]
146      !!-------------------------------------------------------------------------------
147      rho_air_vctr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
148   END FUNCTION rho_air_vctr
149
150   FUNCTION rho_air_sclr( ptak, pqa, pslp )
151      !!-------------------------------------------------------------------------------
152      !!                           ***  FUNCTION rho_air_sclr  ***
153      !!
154      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
155      !!
156      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
157      !!-------------------------------------------------------------------------------
158      REAL(wp), INTENT(in) :: ptak           ! air temperature             [K]
159      REAL(wp), INTENT(in) :: pqa            ! air specific humidity   [kg/kg]
160      REAL(wp), INTENT(in) :: pslp           ! pressure in                [Pa]
161      REAL(wp)             :: rho_air_sclr   ! density of moist air   [kg/m^3]
162      !!-------------------------------------------------------------------------------
163      rho_air_sclr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
164   END FUNCTION rho_air_sclr
165
166
167
168   FUNCTION visc_air(ptak)
169      !!----------------------------------------------------------------------------------
170      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
171      !!
172      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
173      !!----------------------------------------------------------------------------------
174      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
175      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
176      !
177      INTEGER  ::   ji, jj      ! dummy loop indices
178      REAL(wp) ::   ztc, ztc2   ! local scalar
179      !!----------------------------------------------------------------------------------
180      !
181      DO jj = 1, jpj
182         DO ji = 1, jpi
183            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
184            ztc2 = ztc*ztc
185            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
186         END DO
187      END DO
188      !
189   END FUNCTION visc_air
190
191   FUNCTION L_vap_vctr( psst )
192      !!---------------------------------------------------------------------------------
193      !!                           ***  FUNCTION L_vap_vctr  ***
194      !!
195      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
196      !!
197      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
198      !!----------------------------------------------------------------------------------
199      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap_vctr   ! latent heat of vaporization   [J/kg]
200      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
201      !!----------------------------------------------------------------------------------
202      !
203      L_vap_vctr = (  2.501_wp - 0.00237_wp * ( psst(:,:) - rt0)  ) * 1.e6_wp
204      !
205   END FUNCTION L_vap_vctr
206
207   FUNCTION L_vap_sclr( psst )
208      !!---------------------------------------------------------------------------------
209      !!                           ***  FUNCTION L_vap_sclr  ***
210      !!
211      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
212      !!
213      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
214      !!----------------------------------------------------------------------------------
215      REAL(wp)             ::   L_vap_sclr   ! latent heat of vaporization   [J/kg]
216      REAL(wp), INTENT(in) ::   psst         ! water temperature                [K]
217      !!----------------------------------------------------------------------------------
218      !
219      L_vap_sclr = (  2.501_wp - 0.00237_wp * ( psst - rt0)  ) * 1.e6_wp
220      !
221   END FUNCTION L_vap_sclr
222
223
224   FUNCTION cp_air_vctr( pqa )
225      !!-------------------------------------------------------------------------------
226      !!                           ***  FUNCTION cp_air_vctr  ***
227      !!
228      !! ** Purpose : Compute specific heat (Cp) of moist air
229      !!
230      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
231      !!-------------------------------------------------------------------------------
232      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
233      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_vctr   ! specific heat of moist air   [J/K/kg]
234      !!-------------------------------------------------------------------------------
235      cp_air_vctr = rCp_dry + rCp_vap * pqa
236   END FUNCTION cp_air_vctr
237
238   FUNCTION cp_air_sclr( pqa )
239      !!-------------------------------------------------------------------------------
240      !!                           ***  FUNCTION cp_air_sclr  ***
241      !!
242      !! ** Purpose : Compute specific heat (Cp) of moist air
243      !!
244      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
245      !!-------------------------------------------------------------------------------
246      REAL(wp), INTENT(in) :: pqa           ! air specific humidity         [kg/kg]
247      REAL(wp)             :: cp_air_sclr   ! specific heat of moist air   [J/K/kg]
248      !!-------------------------------------------------------------------------------
249      cp_air_sclr = rCp_dry + rCp_vap * pqa
250   END FUNCTION cp_air_sclr
251
252
253
254
255
256   FUNCTION gamma_moist_vctr( ptak, pqa )
257      !!----------------------------------------------------------------------------------
258      !!                           ***  FUNCTION gamma_moist_vctr  ***
259      !!
260      !! ** Purpose : Compute the moist adiabatic lapse-rate.
261      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
262      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
263      !!
264      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
265      !!----------------------------------------------------------------------------------
266      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
267      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
268      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist_vctr   ! moist adiabatic lapse-rate
269      !
270      INTEGER  ::   ji, jj         ! dummy loop indices
271      !!----------------------------------------------------------------------------------
272      DO jj = 1, jpj
273         DO ji = 1, jpi
274            gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
275         END DO
276      END DO
277   END FUNCTION gamma_moist_vctr
278
279   FUNCTION gamma_moist_sclr( ptak, pqa )
280      !!----------------------------------------------------------------------------------
281      !! ** Purpose : Compute the moist adiabatic lapse-rate.
282      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
283      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
284      !!
285      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
286      !!----------------------------------------------------------------------------------
287      REAL(wp)             :: gamma_moist_sclr
288      REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
289      !
290      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
291      !!----------------------------------------------------------------------------------
292      zta = MAX( ptak,  180._wp) ! prevents screw-up over masked regions where field == 0.
293      zqa = MAX( pqa,  1.E-6_wp) !    "                   "                     "
294      !!
295      zwa = zqa / (1._wp - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
296      ziRT = 1._wp / (R_dry*zta)    ! 1/RT
297      gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
298      !!
299   END FUNCTION gamma_moist_sclr
300
301   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
302      !!------------------------------------------------------------------------
303      !!
304      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
305      !!  specific humidity, and frictional scales u*, t* and q*
306      !!
307      !! Author: L. Brodeau, June 2016 / AeroBulk
308      !!         (https://github.com/brodeau/aerobulk/)
309      !!------------------------------------------------------------------------
310      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
311      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
312         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
313         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
314      !
315      INTEGER  ::   ji, jj         ! dummy loop indices
316      REAL(wp) ::     zqa          ! local scalar
317      !!-------------------------------------------------------------------
318      !
319      DO jj = 1, jpj
320         DO ji = 1, jpi
321            !
322            zqa = (1._wp + rctv0*pqa(ji,jj))
323            !
324            ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
325            !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
326            !                      or
327            !  b/  -u* [ theta*              + 0.61 theta q* ]
328            !
329            One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
330               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
331            !
332         END DO
333      END DO
334      !
335      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
336      !
337   END FUNCTION One_on_L
338
339   FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub )
340      !!----------------------------------------------------------------------------------
341      !! Bulk Richardson number according to "wide-spread equation"...
342      !!
343      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
344      !!----------------------------------------------------------------------------------
345      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk
346      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m]
347      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst  ! SST                                   [K]
348      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K]
349      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg]
350      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg]
351      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub   ! bulk wind speed                     [m/s]
352      !
353      INTEGER  ::   ji, jj                                ! dummy loop indices
354      REAL(wp) ::   zqa, zta, zgamma, zdth_v, ztv, zsstv  ! local scalars
355      !!-------------------------------------------------------------------
356      !
357      DO jj = 1, jpj
358         DO ji = 1, jpi
359            !
360            zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer...
361            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
362            zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer
363            zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer
364            !
365            zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!)
366            !
367            zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature"
368            !
369            ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer
370            !
371            Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk
372            !
373         END DO
374      END DO
375   END FUNCTION Ri_bulk
376
377
378   FUNCTION e_sat_vctr(ptak)
379      !!**************************************************
380      !! ptak:     air temperature [K]
381      !! e_sat:  water vapor at saturation [Pa]
382      !!
383      !! Recommended by WMO
384      !!
385      !! Goff, J. A., 1957: Saturation pressure of water on the new kelvin
386      !! temperature scale. Transactions of the American society of heating
387      !! and ventilating engineers, 347–354.
388      !!
389      !! rt0 should be 273.16 (triple point of water) and not 273.15 like here
390      !!**************************************************
391
392      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa]
393      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak    !: temperature (K)
394
395      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp
396
397      ALLOCATE ( ztmp(jpi,jpj) )
398
399      ztmp(:,:) = rtt0/ptak(:,:)
400
401      e_sat_vctr = 100.*( 10.**(10.79574*(1. - ztmp) - 5.028*LOG10(ptak/rtt0)         &
402         &       + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak/rtt0 - 1.)) )   &
403         &       + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
404
405      DEALLOCATE ( ztmp )
406
407   END FUNCTION e_sat_vctr
408
409
410   FUNCTION e_sat_sclr( ptak )
411      !!----------------------------------------------------------------------------------
412      !!                   ***  FUNCTION e_sat_sclr  ***
413      !!                  < SCALAR argument version >
414      !! ** Purpose : water vapor at saturation in [Pa]
415      !!              Based on accurate estimate by Goff, 1957
416      !!
417      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
418      !!
419      !!    Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here
420      !!----------------------------------------------------------------------------------
421      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K]
422      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg]
423      !
424      REAL(wp) ::   zta, ztmp   ! local scalar
425      !!----------------------------------------------------------------------------------
426      !
427      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
428      ztmp = rt0 / zta
429      !
430      ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957)
431      e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        &
432         &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  &
433         &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
434      !
435   END FUNCTION e_sat_sclr
436
437
438   FUNCTION q_sat( ptak, pslp )
439      !!----------------------------------------------------------------------------------
440      !!                           ***  FUNCTION q_sat  ***
441      !!
442      !! ** Purpose : Specific humidity at saturation in [kg/kg]
443      !!              Based on accurate estimate of "e_sat"
444      !!              aka saturation water vapor (Goff, 1957)
445      !!
446      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
447      !!----------------------------------------------------------------------------------
448      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
449      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
450      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
451      !
452      INTEGER  ::   ji, jj         ! dummy loop indices
453      REAL(wp) ::   ze_sat   ! local scalar
454      !!----------------------------------------------------------------------------------
455      !
456      DO jj = 1, jpj
457         DO ji = 1, jpi
458            !
459            ze_sat =  e_sat_sclr( ptak(ji,jj) )
460            !
461            q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat )
462            !
463         END DO
464      END DO
465      !
466   END FUNCTION q_sat
467
468   FUNCTION q_air_rh(prha, ptak, pslp)
469      !!----------------------------------------------------------------------------------
470      !! Specific humidity of air out of Relative Humidity
471      !!
472      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
473      !!----------------------------------------------------------------------------------
474      REAL(wp), DIMENSION(jpi,jpj)             :: q_air_rh
475      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha        !: relative humidity      [fraction, not %!!!]
476      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak        !: air temperature        [K]
477      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp        !: atmospheric pressure   [Pa]
478      !
479      INTEGER  ::   ji, jj      ! dummy loop indices
480      REAL(wp) ::   ze      ! local scalar
481      !!----------------------------------------------------------------------------------
482      !
483      DO jj = 1, jpj
484         DO ji = 1, jpi
485            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
486            q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze)
487         END DO
488      END DO
489      !
490   END FUNCTION q_air_rh
491
492
493   SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, pslp, prlw, &
494      &                         pQns, pTau,  &
495      &                         Qlat)
496      !!----------------------------------------------------------------------------------
497      !! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw"
498      !!          and the module of the wind stress => pTau = Tau
499      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
500      !!----------------------------------------------------------------------------------
501      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
502      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
503      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
504      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
505      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
506      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u*
507      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t*
508      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqst ! q*
509      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
510      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
511      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
512      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2]
513      !
514      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
515      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
516      !
517      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat
518      !
519      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zTs2, zz0, &
520         &        zQlat, zQsen, zQlw
521      INTEGER  ::   ji, jj     ! dummy loop indices
522      !!----------------------------------------------------------------------------------
523      DO jj = 1, jpj
524         DO ji = 1, jpi
525
526            zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
527            zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
528            zz0 = pust(ji,jj)/pUb(ji,jj)
529            zCd = zz0*zz0
530            zCh = zz0*ptst(ji,jj)/zdt
531            zCe = zz0*pqst(ji,jj)/zdq
532
533            CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
534               &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), &
535               &              pTau(ji,jj), zQsen, zQlat )
536
537            zTs2  = pTs(ji,jj)*pTs(ji,jj)
538            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
539
540            pQns(ji,jj) = zQlat + zQsen + zQlw
541
542            IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
543         END DO
544      END DO
545   END SUBROUTINE UPDATE_QNSOL_TAU
546
547
548   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, &
549      &                                 pTau, pQsen, pQlat,  pEvap, prhoa )
550      !!----------------------------------------------------------------------------------
551      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
552      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
553      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
554      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
555      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
556      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCd
557      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCh
558      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCe
559      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
560      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
561      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
562      !!
563      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
564      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen !  [W/m^2]
565      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2]
566      !!
567      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
568      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
569      !!
570      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
571      INTEGER  :: ji, jj, jq     ! dummy loop indices
572      !!----------------------------------------------------------------------------------
573      DO jj = 1, jpj
574         DO ji = 1, jpi
575
576            !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
577            ztaa = pTa(ji,jj) ! first guess...
578            DO jq = 1, 4
579               zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) )
580               ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder...
581            END DO
582            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj))
583            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
584
585            zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10
586
587            pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module
588
589            zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj))
590            pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
591            pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
592
593            IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
594            IF ( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
595
596         END DO
597      END DO
598   END SUBROUTINE BULK_FORMULA_VCTR
599
600
601   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, &
602      &                                 pTau, pQsen, pQlat,  pEvap, prhoa )
603      !!----------------------------------------------------------------------------------
604      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
605      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
606      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
607      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
608      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
609      REAL(wp), INTENT(in)  :: pCd
610      REAL(wp), INTENT(in)  :: pCh
611      REAL(wp), INTENT(in)  :: pCe
612      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
613      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
614      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
615      !!
616      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
617      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
618      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
619      !!
620      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
621      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
622      !!
623      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
624      INTEGER  :: jq
625      !!----------------------------------------------------------------------------------
626
627      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
628      ztaa = pTa ! first guess...
629      DO jq = 1, 4
630         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )
631         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
632      END DO
633      zrho = rho_air(ztaa, pqa, pslp)
634      zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
635
636      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
637
638      pTau = zUrho * pCd * pwnd ! Wind stress module
639
640      zevap = zUrho * pCe * (pqa - pqs)
641      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
642      pQlat = L_vap(pTs) * zevap
643
644      IF ( PRESENT(pEvap) ) pEvap = - zevap
645      IF ( PRESENT(prhoa) ) prhoa = zrho
646
647   END SUBROUTINE BULK_FORMULA_SCLR
648
649
650
651
652   FUNCTION alpha_sw_vctr( psst )
653      !!---------------------------------------------------------------------------------
654      !!                           ***  FUNCTION alpha_sw_vctr  ***
655      !!
656      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
657      !!
658      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
659      !!----------------------------------------------------------------------------------
660      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! thermal expansion coefficient of sea-water [1/K]
661      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
662      !!----------------------------------------------------------------------------------
663      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
664   END FUNCTION alpha_sw_vctr
665
666   FUNCTION alpha_sw_sclr( psst )
667      !!---------------------------------------------------------------------------------
668      !!                           ***  FUNCTION alpha_sw_sclr  ***
669      !!
670      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
671      !!
672      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
673      !!----------------------------------------------------------------------------------
674      REAL(wp)             ::   alpha_sw_sclr   ! thermal expansion coefficient of sea-water [1/K]
675      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K]
676      !!----------------------------------------------------------------------------------
677      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
678   END FUNCTION alpha_sw_sclr
679
680
681
682   !!======================================================================
683END MODULE sbcblk_phy
Note: See TracBrowser for help on using the repository browser.