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_ASINTER-01-05_merged/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_phy.F90 @ 12015

Last change on this file since 12015 was 12015, checked in by gsamson, 4 years ago

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