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

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

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