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/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcblk_phy.F90 @ 12724

Last change on this file since 12724 was 12724, checked in by techene, 4 years ago

branch KERNEL-06 : merge with trunk@12698 #2385 - in duplcated files : changes to comply to the new trunk variables and some loop bug fixes
File size: 35.2 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 longwave 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   !! * Substitutions
103#  include "do_loop_substitute.h90"
104   !!----------------------------------------------------------------------
105   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
106   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
107   !! Software governed by the CeCILL license (see ./LICENSE)
108   !!----------------------------------------------------------------------
109CONTAINS
110
111   FUNCTION virt_temp( pta, pqa )
112      !!------------------------------------------------------------------------
113      !!
114      !! Compute the (absolute/potential) virtual temperature, knowing the
115      !! (absolute/potential) temperature and specific humidity
116      !!
117      !! If input temperature is absolute then output vitual temperature is absolute
118      !! If input temperature is potential then output vitual temperature is potential
119      !!
120      !! Author: L. Brodeau, June 2019 / AeroBulk
121      !!         (https://github.com/brodeau/aerobulk/)
122      !!------------------------------------------------------------------------
123      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp         !: 1./(Monin Obukhov length) [m^-1]
124      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta,  &  !: absolute or potetntial air temperature [K]
125         &                                        pqa      !: specific humidity of air   [kg/kg]
126      !!-------------------------------------------------------------------
127      !
128      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
129      !!
130      !! This is exactly the same sing that:
131      !! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa))
132      !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
133      !
134   END FUNCTION virt_temp
135
136   FUNCTION rho_air_vctr( ptak, pqa, pslp )
137      !!-------------------------------------------------------------------------------
138      !!                           ***  FUNCTION rho_air_vctr  ***
139      !!
140      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
141      !!
142      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
143      !!-------------------------------------------------------------------------------
144      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
145      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
146      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
147      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air_vctr   ! density of moist air   [kg/m^3]
148      !!-------------------------------------------------------------------------------
149      rho_air_vctr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
150   END FUNCTION rho_air_vctr
151
152   FUNCTION rho_air_sclr( ptak, pqa, pslp )
153      !!-------------------------------------------------------------------------------
154      !!                           ***  FUNCTION rho_air_sclr  ***
155      !!
156      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
157      !!
158      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
159      !!-------------------------------------------------------------------------------
160      REAL(wp), INTENT(in) :: ptak           ! air temperature             [K]
161      REAL(wp), INTENT(in) :: pqa            ! air specific humidity   [kg/kg]
162      REAL(wp), INTENT(in) :: pslp           ! pressure in                [Pa]
163      REAL(wp)             :: rho_air_sclr   ! density of moist air   [kg/m^3]
164      !!-------------------------------------------------------------------------------
165      rho_air_sclr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
166   END FUNCTION rho_air_sclr
167
168
169
170   FUNCTION visc_air(ptak)
171      !!----------------------------------------------------------------------------------
172      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
173      !!
174      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
175      !!----------------------------------------------------------------------------------
176      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
177      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
178      !
179      INTEGER  ::   ji, jj      ! dummy loop indices
180      REAL(wp) ::   ztc, ztc2   ! local scalar
181      !!----------------------------------------------------------------------------------
182      !
183      DO_2D_11_11
184         ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
185         ztc2 = ztc*ztc
186         visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
187      END_2D
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_2D_11_11
273         gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
274      END_2D
275   END FUNCTION gamma_moist_vctr
276
277   FUNCTION gamma_moist_sclr( ptak, pqa )
278      !!----------------------------------------------------------------------------------
279      !! ** Purpose : Compute the moist adiabatic lapse-rate.
280      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
281      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
282      !!
283      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
284      !!----------------------------------------------------------------------------------
285      REAL(wp)             :: gamma_moist_sclr
286      REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
287      !
288      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
289      !!----------------------------------------------------------------------------------
290      zta = MAX( ptak,  180._wp) ! prevents screw-up over masked regions where field == 0.
291      zqa = MAX( pqa,  1.E-6_wp) !    "                   "                     "
292      !!
293      zwa = zqa / (1._wp - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
294      ziRT = 1._wp / (R_dry*zta)    ! 1/RT
295      gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
296      !!
297   END FUNCTION gamma_moist_sclr
298
299   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
300      !!------------------------------------------------------------------------
301      !!
302      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
303      !!  specific humidity, and frictional scales u*, t* and q*
304      !!
305      !! Author: L. Brodeau, June 2016 / AeroBulk
306      !!         (https://github.com/brodeau/aerobulk/)
307      !!------------------------------------------------------------------------
308      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
309      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
310         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
311         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
312      !
313      INTEGER  ::   ji, jj         ! dummy loop indices
314      REAL(wp) ::     zqa          ! local scalar
315      !!-------------------------------------------------------------------
316      !
317      DO_2D_11_11
318         !
319         zqa = (1._wp + rctv0*pqa(ji,jj))
320         !
321         ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
322         !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
323         !                      or
324         !  b/  -u* [ theta*              + 0.61 theta q* ]
325         !
326         One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
327            &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
328         !
329      END_2D
330      !
331      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
332      !
333   END FUNCTION One_on_L
334
335   FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub )
336      !!----------------------------------------------------------------------------------
337      !! Bulk Richardson number according to "wide-spread equation"...
338      !!
339      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
340      !!----------------------------------------------------------------------------------
341      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk
342      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m]
343      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst  ! SST                                   [K]
344      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K]
345      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg]
346      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg]
347      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub   ! bulk wind speed                     [m/s]
348      !
349      INTEGER  ::   ji, jj                                ! dummy loop indices
350      REAL(wp) ::   zqa, zta, zgamma, zdth_v, ztv, zsstv  ! local scalars
351      !!-------------------------------------------------------------------
352      !
353      DO_2D_11_11
354         !
355         zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer...
356         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
357         zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer
358         zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer
359         !
360         zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!)
361         !
362         zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature"
363         !
364         ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer
365         !
366         Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk
367         !
368      END_2D
369   END FUNCTION Ri_bulk
370
371
372   FUNCTION e_sat_vctr(ptak)
373      !!**************************************************
374      !! ptak:     air temperature [K]
375      !! e_sat:  water vapor at saturation [Pa]
376      !!
377      !! Recommended by WMO
378      !!
379      !! Goff, J. A., 1957: Saturation pressure of water on the new kelvin
380      !! temperature scale. Transactions of the American society of heating
381      !! and ventilating engineers, 347–354.
382      !!
383      !! rt0 should be 273.16 (triple point of water) and not 273.15 like here
384      !!**************************************************
385
386      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa]
387      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak    !: temperature (K)
388
389      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp
390
391      ALLOCATE ( ztmp(jpi,jpj) )
392
393      ztmp(:,:) = rtt0/ptak(:,:)
394
395      e_sat_vctr = 100.*( 10.**(10.79574*(1. - ztmp) - 5.028*LOG10(ptak/rtt0)         &
396         &       + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak/rtt0 - 1.)) )   &
397         &       + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
398
399      DEALLOCATE ( ztmp )
400
401   END FUNCTION e_sat_vctr
402
403
404   FUNCTION e_sat_sclr( ptak )
405      !!----------------------------------------------------------------------------------
406      !!                   ***  FUNCTION e_sat_sclr  ***
407      !!                  < SCALAR argument version >
408      !! ** Purpose : water vapor at saturation in [Pa]
409      !!              Based on accurate estimate by Goff, 1957
410      !!
411      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
412      !!
413      !!    Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here
414      !!----------------------------------------------------------------------------------
415      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K]
416      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg]
417      !
418      REAL(wp) ::   zta, ztmp   ! local scalar
419      !!----------------------------------------------------------------------------------
420      !
421      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
422      ztmp = rt0 / zta
423      !
424      ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957)
425      e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        &
426         &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  &
427         &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
428      !
429   END FUNCTION e_sat_sclr
430
431
432   FUNCTION q_sat( ptak, pslp )
433      !!----------------------------------------------------------------------------------
434      !!                           ***  FUNCTION q_sat  ***
435      !!
436      !! ** Purpose : Specific humidity at saturation in [kg/kg]
437      !!              Based on accurate estimate of "e_sat"
438      !!              aka saturation water vapor (Goff, 1957)
439      !!
440      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
441      !!----------------------------------------------------------------------------------
442      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
443      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
444      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
445      !
446      INTEGER  ::   ji, jj         ! dummy loop indices
447      REAL(wp) ::   ze_sat   ! local scalar
448      !!----------------------------------------------------------------------------------
449      !
450      DO_2D_11_11
451         !
452         ze_sat =  e_sat_sclr( ptak(ji,jj) )
453         !
454         q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat )
455         !
456      END_2D
457      !
458   END FUNCTION q_sat
459
460   FUNCTION q_air_rh(prha, ptak, pslp)
461      !!----------------------------------------------------------------------------------
462      !! Specific humidity of air out of Relative Humidity
463      !!
464      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
465      !!----------------------------------------------------------------------------------
466      REAL(wp), DIMENSION(jpi,jpj)             :: q_air_rh
467      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha        !: relative humidity      [fraction, not %!!!]
468      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak        !: air temperature        [K]
469      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp        !: atmospheric pressure   [Pa]
470      !
471      INTEGER  ::   ji, jj      ! dummy loop indices
472      REAL(wp) ::   ze      ! local scalar
473      !!----------------------------------------------------------------------------------
474      !
475      DO_2D_11_11
476         ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
477         q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze)
478      END_2D
479      !
480   END FUNCTION q_air_rh
481
482
483   SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, pslp, prlw, &
484      &                         pQns, pTau,  &
485      &                         Qlat)
486      !!----------------------------------------------------------------------------------
487      !! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw"
488      !!          and the module of the wind stress => pTau = Tau
489      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
490      !!----------------------------------------------------------------------------------
491      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
492      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
493      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
494      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
495      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
496      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u*
497      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t*
498      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqst ! q*
499      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
500      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
501      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
502      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2]
503      !
504      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
505      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
506      !
507      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat
508      !
509      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zTs2, zz0, &
510         &        zQlat, zQsen, zQlw
511      INTEGER  ::   ji, jj     ! dummy loop indices
512      !!----------------------------------------------------------------------------------
513      DO_2D_11_11
514
515         zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
516         zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
517         zz0 = pust(ji,jj)/pUb(ji,jj)
518         zCd = zz0*zz0
519         zCh = zz0*ptst(ji,jj)/zdt
520         zCe = zz0*pqst(ji,jj)/zdq
521
522         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
523            &                    zCd, zCh, zCe,                                       &
524            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                &
525            &                    pTau(ji,jj), zQsen, zQlat )
526         
527         zTs2  = pTs(ji,jj)*pTs(ji,jj)
528         zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
529
530         pQns(ji,jj) = zQlat + zQsen + zQlw
531
532         IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
533      END_2D
534   END SUBROUTINE UPDATE_QNSOL_TAU
535
536
537   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
538      &                          pCd, pCh, pCe,           &
539      &                          pwnd, pUb, pslp,         &
540      &                          pTau, pQsen, pQlat,      &
541      &                          pEvap, prhoa, pfact_evap )
542      !!----------------------------------------------------------------------------------
543      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
544      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
545      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
546      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
547      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
548      REAL(wp), INTENT(in)  :: pCd
549      REAL(wp), INTENT(in)  :: pCh
550      REAL(wp), INTENT(in)  :: pCe
551      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
552      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
553      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
554      !!
555      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
556      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
557      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
558      !!
559      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
560      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
561      REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
562      !!
563      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
564      INTEGER  :: jq
565      !!----------------------------------------------------------------------------------
566      zfact_evap = 1._wp
567      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
568     
569      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
570      ztaa = pTa ! first guess...
571      DO jq = 1, 4
572         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )  !LOLO: why not "0.5*(pqs+pqa)" rather then "pqa" ???
573         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
574      END DO
575      zrho = rho_air(ztaa, pqa, pslp)
576      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!
577
578      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
579
580      pTau = zUrho * pCd * pwnd ! Wind stress module
581
582      zevap = zUrho * pCe * (pqa - pqs)
583      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
584      pQlat = L_vap(pTs) * zevap
585
586      IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap
587      IF( PRESENT(prhoa) ) prhoa = zrho
588
589   END SUBROUTINE BULK_FORMULA_SCLR
590
591   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
592      &                          pCd, pCh, pCe,           &
593      &                          pwnd, pUb, pslp,         &
594      &                          pTau, pQsen, pQlat,      & 
595      &                          pEvap, prhoa, pfact_evap )     
596      !!----------------------------------------------------------------------------------
597      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
598      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
599      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
600      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
601      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
602      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCd
603      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCh
604      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCe
605      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
606      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
607      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
608      !!
609      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
610      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen !  [W/m^2]
611      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2]
612      !!
613      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
614      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
615      REAL(wp),                     INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
616      !!
617      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
618      INTEGER  :: ji, jj
619      !!----------------------------------------------------------------------------------
620      zfact_evap = 1._wp
621      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
622
623      DO_2D_11_11
624
625         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
626            &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  &
627            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                &
628            &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             &
629            &                    pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap       )
630
631         IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
632         IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
633   
634      END_2D
635   END SUBROUTINE BULK_FORMULA_VCTR
636
637
638   FUNCTION alpha_sw_vctr( psst )
639      !!---------------------------------------------------------------------------------
640      !!                           ***  FUNCTION alpha_sw_vctr  ***
641      !!
642      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
643      !!
644      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
645      !!----------------------------------------------------------------------------------
646      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! thermal expansion coefficient of sea-water [1/K]
647      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
648      !!----------------------------------------------------------------------------------
649      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
650   END FUNCTION alpha_sw_vctr
651
652   FUNCTION alpha_sw_sclr( psst )
653      !!---------------------------------------------------------------------------------
654      !!                           ***  FUNCTION alpha_sw_sclr  ***
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)             ::   alpha_sw_sclr   ! thermal expansion coefficient of sea-water [1/K]
661      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K]
662      !!----------------------------------------------------------------------------------
663      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
664   END FUNCTION alpha_sw_sclr
665
666
667
668   !!======================================================================
669END MODULE sbcblk_phy
Note: See TracBrowser for help on using the repository browser.