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 @ 12162

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

Catching up with "dev_r12072_MERGE_OPTION2_2019"

File size: 35.5 KB
Line 
1MODULE sbcblk_phy
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_phy  ***
4   !! A set of functions to compute air themodynamics parameters
5   !!                     needed by Aerodynamic Bulk Formulas
6   !!=====================================================================
7   !!            4.0  !  2019 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
8   !!----------------------------------------------------------------------
9
10   !!   virt_temp     : virtual (aka sensible) temperature (potential or absolute)
11   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
12   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature
13   !!   L_vap         : latent heat of vaporization of water as a function of temperature
14   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
15   !!   gamma_moist   : adiabatic lapse-rate of moist air
16   !!   One_on_L      : 1. / ( Monin-Obukhov length )
17   !!   Ri_bulk       : bulk Richardson number aka BRN
18   !!   q_sat         : saturation humidity as a function of SLP and temperature
19   !!   q_air_rh      : specific humidity as a function of RH (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   !!----------------------------------------------------------------------
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), &
534               &              zCd, zCh, zCe,                                        &
535               &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                 &
536               &              pTau(ji,jj), zQsen, zQlat )
537
538            zTs2  = pTs(ji,jj)*pTs(ji,jj)
539            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
540
541            pQns(ji,jj) = zQlat + zQsen + zQlw
542
543            IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
544         END DO
545      END DO
546   END SUBROUTINE UPDATE_QNSOL_TAU
547
548
549   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa,  &
550      &                          pCd, pCh, pCe,            &
551      &                          pwnd, pUb, pslp,          &
552      &                          pTau, pQsen, pQlat,  pEvap, prhoa )
553      !!----------------------------------------------------------------------------------
554      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
555      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
556      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
557      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
558      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
559      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCd
560      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCh
561      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCe
562      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
563      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
564      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
565      !!
566      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
567      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen !  [W/m^2]
568      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2]
569      !!
570      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
571      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
572      !!
573      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
574      INTEGER  :: ji, jj, jq     ! dummy loop indices
575      !!----------------------------------------------------------------------------------
576      DO jj = 1, jpj
577         DO ji = 1, jpi
578
579            !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
580            ztaa = pTa(ji,jj) ! first guess...
581            DO jq = 1, 4
582               zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) )
583               ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder...
584            END DO
585            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj))
586            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!
587
588            zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10
589
590            pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module
591
592            zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj))
593            pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
594            pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
595
596            IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
597            IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
598
599         END DO
600      END DO
601   END SUBROUTINE BULK_FORMULA_VCTR
602
603
604   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
605      &                          pCd, pCh, pCe,           &
606      &                          pwnd, pUb, pslp,         &
607      &                          pTau, pQsen, pQlat,  pEvap, prhoa )
608      !!----------------------------------------------------------------------------------
609      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
610      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
611      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
612      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
613      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
614      REAL(wp), INTENT(in)  :: pCd
615      REAL(wp), INTENT(in)  :: pCh
616      REAL(wp), INTENT(in)  :: pCe
617      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
618      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
619      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
620      !!
621      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
622      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
623      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
624      !!
625      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
626      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
627      !!
628      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
629      INTEGER  :: jq
630      !!----------------------------------------------------------------------------------
631
632      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
633      ztaa = pTa ! first guess...
634      DO jq = 1, 4
635         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )
636         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
637      END DO
638      zrho = rho_air(ztaa, pqa, pslp)
639      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!
640
641      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
642
643      pTau = zUrho * pCd * pwnd ! Wind stress module
644
645      zevap = zUrho * pCe * (pqa - pqs)
646      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
647      pQlat = L_vap(pTs) * zevap
648
649      IF( PRESENT(pEvap) ) pEvap = - zevap
650      IF( PRESENT(prhoa) ) prhoa = zrho
651
652   END SUBROUTINE BULK_FORMULA_SCLR
653
654
655
656
657   FUNCTION alpha_sw_vctr( psst )
658      !!---------------------------------------------------------------------------------
659      !!                           ***  FUNCTION alpha_sw_vctr  ***
660      !!
661      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
662      !!
663      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
664      !!----------------------------------------------------------------------------------
665      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! thermal expansion coefficient of sea-water [1/K]
666      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
667      !!----------------------------------------------------------------------------------
668      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
669   END FUNCTION alpha_sw_vctr
670
671   FUNCTION alpha_sw_sclr( psst )
672      !!---------------------------------------------------------------------------------
673      !!                           ***  FUNCTION alpha_sw_sclr  ***
674      !!
675      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
676      !!
677      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
678      !!----------------------------------------------------------------------------------
679      REAL(wp)             ::   alpha_sw_sclr   ! thermal expansion coefficient of sea-water [1/K]
680      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K]
681      !!----------------------------------------------------------------------------------
682      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
683   END FUNCTION alpha_sw_sclr
684
685
686
687   !!======================================================================
688END MODULE sbcblk_phy
Note: See TracBrowser for help on using the repository browser.