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

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

Last change on this file since 11841 was 11841, checked in by laurent, 4 years ago

Better use of function "bulk_formula()" (in sbcblk_phy.F90)
File size: 35.7 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      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
272      !!----------------------------------------------------------------------------------
273      !
274      DO jj = 1, jpj
275         DO ji = 1, jpi
276            zta = MAX( ptak(ji,jj),  180._wp) ! prevents screw-up over masked regions where field == 0.
277            zqa = MAX( pqa(ji,jj),  1.E-6_wp) !    "                   "                     "
278            !
279            zwa = zqa / (1. - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
280            ziRT = 1._wp/(R_dry*zta)    ! 1/RT
281            gamma_moist_vctr(ji,jj) = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
282         END DO
283      END DO
284      !
285   END FUNCTION gamma_moist_vctr
286
287   FUNCTION gamma_moist_sclr( ptak, pqa )
288      !!----------------------------------------------------------------------------------
289      !! ** Purpose : Compute the moist adiabatic lapse-rate.
290      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
291      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
292      !!
293      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
294      !!----------------------------------------------------------------------------------
295      REAL(wp)             :: gamma_moist_sclr
296      REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
297      !
298      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
299      !!----------------------------------------------------------------------------------
300      zta = MAX( ptak,  180._wp) ! prevents screw-up over masked regions where field == 0.
301      zqa = MAX( pqa,  1.E-6_wp) !    "                   "                     "
302      !!
303      zwa = zqa / (1._wp - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
304      ziRT = 1._wp / (R_dry*zta)    ! 1/RT
305      gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
306      !!
307   END FUNCTION gamma_moist_sclr
308
309   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
310      !!------------------------------------------------------------------------
311      !!
312      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
313      !!  specific humidity, and frictional scales u*, t* and q*
314      !!
315      !! Author: L. Brodeau, June 2016 / AeroBulk
316      !!         (https://github.com/brodeau/aerobulk/)
317      !!------------------------------------------------------------------------
318      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
319      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
320         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
321         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
322      !
323      INTEGER  ::   ji, jj         ! dummy loop indices
324      REAL(wp) ::     zqa          ! local scalar
325      !!-------------------------------------------------------------------
326      !
327      DO jj = 1, jpj
328         DO ji = 1, jpi
329            !
330            zqa = (1._wp + rctv0*pqa(ji,jj))
331            !
332            ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
333            !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
334            !                      or
335            !  b/  -u* [ theta*              + 0.61 theta q* ]
336            !
337            One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
338               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
339            !
340         END DO
341      END DO
342      !
343      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
344      !
345   END FUNCTION One_on_L
346
347   FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub )
348      !!----------------------------------------------------------------------------------
349      !! Bulk Richardson number according to "wide-spread equation"...
350      !!
351      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
352      !!----------------------------------------------------------------------------------
353      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk
354      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m]
355      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst  ! SST                                   [K]
356      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K]
357      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg]
358      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg]
359      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub   ! bulk wind speed                     [m/s]
360      !
361      INTEGER  ::   ji, jj                                ! dummy loop indices
362      REAL(wp) ::   zqa, zta, zgamma, zdth_v, ztv, zsstv  ! local scalars
363      !!-------------------------------------------------------------------
364      !
365      DO jj = 1, jpj
366         DO ji = 1, jpi
367            !
368            zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer...
369            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
370            zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer
371            zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer
372            !
373            zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!)
374            !
375            zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature"
376            !
377            ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer
378            !
379            Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk
380            !
381         END DO
382      END DO
383   END FUNCTION Ri_bulk
384
385
386   FUNCTION e_sat_vctr(ptak)
387      !!**************************************************
388      !! ptak:     air temperature [K]
389      !! e_sat:  water vapor at saturation [Pa]
390      !!
391      !! Recommended by WMO
392      !!
393      !! Goff, J. A., 1957: Saturation pressure of water on the new kelvin
394      !! temperature scale. Transactions of the American society of heating
395      !! and ventilating engineers, 347–354.
396      !!
397      !! rt0 should be 273.16 (triple point of water) and not 273.15 like here
398      !!**************************************************
399
400      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa]
401      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak    !: temperature (K)
402
403      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp
404
405      ALLOCATE ( ztmp(jpi,jpj) )
406
407      ztmp(:,:) = rtt0/ptak(:,:)
408
409      e_sat_vctr = 100.*( 10.**(10.79574*(1. - ztmp) - 5.028*LOG10(ptak/rtt0)         &
410         &       + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak/rtt0 - 1.)) )   &
411         &       + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
412
413      DEALLOCATE ( ztmp )
414
415   END FUNCTION e_sat_vctr
416
417
418   FUNCTION e_sat_sclr( ptak )
419      !!----------------------------------------------------------------------------------
420      !!                   ***  FUNCTION e_sat_sclr  ***
421      !!                  < SCALAR argument version >
422      !! ** Purpose : water vapor at saturation in [Pa]
423      !!              Based on accurate estimate by Goff, 1957
424      !!
425      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
426      !!
427      !!    Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here
428      !!----------------------------------------------------------------------------------
429      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K]
430      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg]
431      !
432      REAL(wp) ::   zta, ztmp   ! local scalar
433      !!----------------------------------------------------------------------------------
434      !
435      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
436      ztmp = rt0 / zta
437      !
438      ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957)
439      e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        &
440         &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  &
441         &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
442      !
443   END FUNCTION e_sat_sclr
444
445
446   FUNCTION q_sat( ptak, pslp )
447      !!----------------------------------------------------------------------------------
448      !!                           ***  FUNCTION q_sat  ***
449      !!
450      !! ** Purpose : Specific humidity at saturation in [kg/kg]
451      !!              Based on accurate estimate of "e_sat"
452      !!              aka saturation water vapor (Goff, 1957)
453      !!
454      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
455      !!----------------------------------------------------------------------------------
456      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
457      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
458      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
459      !
460      INTEGER  ::   ji, jj         ! dummy loop indices
461      REAL(wp) ::   ze_sat   ! local scalar
462      !!----------------------------------------------------------------------------------
463      !
464      DO jj = 1, jpj
465         DO ji = 1, jpi
466            !
467            ze_sat =  e_sat_sclr( ptak(ji,jj) )
468            !
469            q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat )
470            !
471         END DO
472      END DO
473      !
474   END FUNCTION q_sat
475
476   FUNCTION q_air_rh(prha, ptak, pslp)
477      !!----------------------------------------------------------------------------------
478      !! Specific humidity of air out of Relative Humidity
479      !!
480      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
481      !!----------------------------------------------------------------------------------
482      REAL(wp), DIMENSION(jpi,jpj)             :: q_air_rh
483      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha        !: relative humidity      [fraction, not %!!!]
484      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak        !: air temperature        [K]
485      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp        !: atmospheric pressure   [Pa]
486      !
487      INTEGER  ::   ji, jj      ! dummy loop indices
488      REAL(wp) ::   ze      ! local scalar
489      !!----------------------------------------------------------------------------------
490      !
491      DO jj = 1, jpj
492         DO ji = 1, jpi
493            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
494            q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze)
495         END DO
496      END DO
497      !
498   END FUNCTION q_air_rh
499
500
501   SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, pslp, prlw, &
502      &                         pQns, pTau,  &
503      &                         Qlat)
504      !!----------------------------------------------------------------------------------
505      !! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw"
506      !!          and the module of the wind stress => pTau = Tau
507      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
508      !!----------------------------------------------------------------------------------
509      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
510      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
511      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
512      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
513      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
514      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u*
515      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t*
516      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqst ! q*
517      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
518      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
519      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
520      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2]
521      !
522      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
523      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
524      !
525      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat
526      !
527      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zTs2, zz0, &
528         &        zQlat, zQsen, zQlw
529      INTEGER  ::   ji, jj     ! dummy loop indices
530      !!----------------------------------------------------------------------------------
531      DO jj = 1, jpj
532         DO ji = 1, jpi
533
534            zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
535            zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
536            zz0 = pust(ji,jj)/pUb(ji,jj)
537            zCd = zz0*zz0
538            zCh = zz0*ptst(ji,jj)/zdt
539            zCe = zz0*pqst(ji,jj)/zdq
540
541            CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
542               &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), &
543               &              pTau(ji,jj), zQsen, zQlat )
544
545            zTs2  = pTs(ji,jj)*pTs(ji,jj)
546            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
547
548            pQns(ji,jj) = zQlat + zQsen + zQlw
549
550            IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
551         END DO
552      END DO
553   END SUBROUTINE UPDATE_QNSOL_TAU
554
555
556   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, &
557      &                                 pTau, pQsen, pQlat,  pEvap, prhoa )
558      !!----------------------------------------------------------------------------------
559      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
560      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
561      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
562      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
563      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
564      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCd
565      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCh
566      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCe
567      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
568      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
569      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
570      !!
571      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
572      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen !  [W/m^2]
573      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2]
574      !!
575      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
576      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
577      !!
578      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
579      INTEGER  :: ji, jj, jq     ! dummy loop indices
580      !!----------------------------------------------------------------------------------
581      DO jj = 1, jpj
582         DO ji = 1, jpi
583
584            !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
585            ztaa = pTa(ji,jj) ! first guess...
586            DO jq = 1, 4
587               zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) )
588               ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder...
589            END DO
590            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj))
591            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!
592
593            zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10
594
595            pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module
596
597            zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj))
598            pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
599            pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
600
601            IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
602            IF ( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
603
604         END DO
605      END DO
606   END SUBROUTINE BULK_FORMULA_VCTR
607
608
609   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, &
610      &                                 pTau, pQsen, pQlat,  pEvap, prhoa )
611      !!----------------------------------------------------------------------------------
612      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
613      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
614      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
615      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
616      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
617      REAL(wp), INTENT(in)  :: pCd
618      REAL(wp), INTENT(in)  :: pCh
619      REAL(wp), INTENT(in)  :: pCe
620      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
621      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
622      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
623      !!
624      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
625      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
626      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
627      !!
628      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
629      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
630      !!
631      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
632      INTEGER  :: jq
633      !!----------------------------------------------------------------------------------
634
635      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
636      ztaa = pTa ! first guess...
637      DO jq = 1, 4
638         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )
639         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
640      END DO
641      zrho = rho_air(ztaa, pqa, pslp)
642      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!
643
644      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
645
646      pTau = zUrho * pCd * pwnd ! Wind stress module
647
648      zevap = zUrho * pCe * (pqa - pqs)
649      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
650      pQlat = L_vap(pTs) * zevap
651
652      IF ( PRESENT(pEvap) ) pEvap = - zevap
653      IF ( PRESENT(prhoa) ) prhoa = zrho
654
655   END SUBROUTINE BULK_FORMULA_SCLR
656
657
658
659
660   FUNCTION alpha_sw_vctr( psst )
661      !!---------------------------------------------------------------------------------
662      !!                           ***  FUNCTION alpha_sw_vctr  ***
663      !!
664      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
665      !!
666      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
667      !!----------------------------------------------------------------------------------
668      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! thermal expansion coefficient of sea-water [1/K]
669      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
670      !!----------------------------------------------------------------------------------
671      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
672   END FUNCTION alpha_sw_vctr
673
674   FUNCTION alpha_sw_sclr( psst )
675      !!---------------------------------------------------------------------------------
676      !!                           ***  FUNCTION alpha_sw_sclr  ***
677      !!
678      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
679      !!
680      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
681      !!----------------------------------------------------------------------------------
682      REAL(wp)             ::   alpha_sw_sclr   ! thermal expansion coefficient of sea-water [1/K]
683      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K]
684      !!----------------------------------------------------------------------------------
685      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
686   END FUNCTION alpha_sw_sclr
687
688
689
690   !!======================================================================
691END MODULE sbcblk_phy
Note: See TracBrowser for help on using the repository browser.