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

Last change on this file since 11615 was 11615, checked in by laurent, 5 years ago

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