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

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

LB: CSWL param now uses NEMO time step (rdt) and previous-t value of t_skin as first guess.
File size: 16.2 KB
Line 
1MODULE sbcblk_phy
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_phy  ***
4   !! A set of functions to compute air themodynamics parameters
5   !!                     needed by Aerodynamic Bulk Formulas
6   !!=====================================================================
7   !!            4.0  !  2019 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
8   !!----------------------------------------------------------------------
9
10   !!   virt_temp     : virtual (aka sensible) temperature (potential or absolute)
11   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
12   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature
13   !!   L_vap         : latent heat of vaporization of water as a function of temperature
14   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
15   !!   gamma_moist   : adiabatic lapse-rate of moist air
16   !!   One_on_L      : 1. / ( Monin-Obukhov length )
17   !!   Ri_bulk       : bulk Richardson number aka BRN
18   !!   q_sat         : saturation humidity as a function of SLP and temperature
19
20   USE dom_oce        ! ocean space and time domain
21   USE phycst         ! physical constants
22
23   IMPLICIT NONE
24   PRIVATE
25
26   INTERFACE gamma_moist
27      MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr
28   END INTERFACE gamma_moist
29
30   PUBLIC virt_temp
31   PUBLIC rho_air
32   PUBLIC visc_air
33   PUBLIC L_vap
34   PUBLIC cp_air
35   PUBLIC gamma_moist
36   PUBLIC One_on_L
37   PUBLIC Ri_bulk
38   PUBLIC q_sat
39
40   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   FUNCTION virt_temp( pta, pqa )
48      !!------------------------------------------------------------------------
49      !!
50      !! Compute the (absolute/potential) virtual temperature, knowing the
51      !! (absolute/potential) temperature and specific humidity
52      !!
53      !! If input temperature is absolute then output vitual temperature is absolute
54      !! If input temperature is potential then output vitual temperature is potential
55      !!
56      !! Author: L. Brodeau, June 2019 / AeroBulk
57      !!         (https://github.com/brodeau/aerobulk/)
58      !!------------------------------------------------------------------------
59      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp         !: 1./(Monin Obukhov length) [m^-1]
60      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta,  &  !: absolute or potetntial air temperature [K]
61         &                                        pqa      !: specific humidity of air   [kg/kg]
62      !!-------------------------------------------------------------------
63      !
64      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
65      !!
66      !! This is exactly the same sing that:
67      !! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa))
68      !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
69      !
70   END FUNCTION virt_temp
71
72   FUNCTION rho_air( ptak, pqa, pslp )
73      !!-------------------------------------------------------------------------------
74      !!                           ***  FUNCTION rho_air  ***
75      !!
76      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
77      !!
78      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
79      !!-------------------------------------------------------------------------------
80      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
81      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
82      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
83      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
84      !!-------------------------------------------------------------------------------
85      !
86      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
87      !
88   END FUNCTION rho_air
89
90   FUNCTION visc_air(ptak)
91      !!----------------------------------------------------------------------------------
92      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
93      !!
94      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
95      !!----------------------------------------------------------------------------------
96      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
97      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
98      !
99      INTEGER  ::   ji, jj      ! dummy loop indices
100      REAL(wp) ::   ztc, ztc2   ! local scalar
101      !!----------------------------------------------------------------------------------
102      !
103      DO jj = 1, jpj
104         DO ji = 1, jpi
105            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
106            ztc2 = ztc*ztc
107            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
108         END DO
109      END DO
110      !
111   END FUNCTION visc_air
112
113   FUNCTION L_vap( psst )
114      !!---------------------------------------------------------------------------------
115      !!                           ***  FUNCTION L_vap  ***
116      !!
117      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
118      !!
119      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
120      !!----------------------------------------------------------------------------------
121      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
122      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
123      !!----------------------------------------------------------------------------------
124      !
125      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
126      !
127   END FUNCTION L_vap
128
129   FUNCTION cp_air( pqa )
130      !!-------------------------------------------------------------------------------
131      !!                           ***  FUNCTION cp_air  ***
132      !!
133      !! ** Purpose : Compute specific heat (Cp) of moist air
134      !!
135      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
136      !!-------------------------------------------------------------------------------
137      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
138      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
139      !!-------------------------------------------------------------------------------
140      !
141      cp_air = rCp_dry + rCp_vap * pqa
142      !
143   END FUNCTION cp_air
144
145   FUNCTION gamma_moist_vctr( ptak, pqa )
146      !!----------------------------------------------------------------------------------
147      !!                           ***  FUNCTION gamma_moist_vctr  ***
148      !!
149      !! ** Purpose : Compute the moist adiabatic lapse-rate.
150      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
151      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
152      !!
153      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
154      !!----------------------------------------------------------------------------------
155      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
156      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
157      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist_vctr   ! moist adiabatic lapse-rate
158      !
159      INTEGER  ::   ji, jj         ! dummy loop indices
160      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
161      !!----------------------------------------------------------------------------------
162      !
163      DO jj = 1, jpj
164         DO ji = 1, jpi
165            zta = MAX( ptak(ji,jj),  180._wp) ! prevents screw-up over masked regions where field == 0.
166            zqa = MAX( pqa(ji,jj),  1.E-6_wp) !    "                   "                     "
167            !
168            zwa = zqa / (1. - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
169            ziRT = 1._wp/(R_dry*zta)    ! 1/RT
170            gamma_moist_vctr(ji,jj) = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
171         END DO
172      END DO
173      !
174   END FUNCTION gamma_moist_vctr
175
176   FUNCTION gamma_moist_sclr( ptak, pqa )
177      !!----------------------------------------------------------------------------------
178      !! ** Purpose : Compute the moist adiabatic lapse-rate.
179      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
180      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
181      !!
182      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
183      !!----------------------------------------------------------------------------------
184      REAL(wp)             :: gamma_moist_sclr
185      REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
186      !
187      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
188      !!----------------------------------------------------------------------------------
189      zta = MAX( ptak,  180._wp) ! prevents screw-up over masked regions where field == 0.
190      zqa = MAX( pqa,  1.E-6_wp) !    "                   "                     "
191      !!
192      zwa = zqa / (1._wp - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
193      ziRT = 1._wp / (R_dry*zta)    ! 1/RT
194      gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
195      !!
196   END FUNCTION gamma_moist_sclr
197
198   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
199      !!------------------------------------------------------------------------
200      !!
201      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
202      !!  specific humidity, and frictional scales u*, t* and q*
203      !!
204      !! Author: L. Brodeau, June 2016 / AeroBulk
205      !!         (https://github.com/brodeau/aerobulk/)
206      !!------------------------------------------------------------------------
207      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
208      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
209         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
210         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
211      !
212      INTEGER  ::   ji, jj         ! dummy loop indices
213      REAL(wp) ::     zqa          ! local scalar
214      !!-------------------------------------------------------------------
215      !
216      DO jj = 1, jpj
217         DO ji = 1, jpi
218            !
219            zqa = (1._wp + rctv0*pqa(ji,jj))
220            !
221            ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
222            !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
223            !                      or
224            !  b/  -u* [ theta*              + 0.61 theta q* ]
225            !
226            One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
227               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
228            !
229         END DO
230      END DO
231      !
232      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
233      !
234   END FUNCTION One_on_L
235
236   FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub )
237      !!----------------------------------------------------------------------------------
238      !! Bulk Richardson number according to "wide-spread equation"...
239      !!
240      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
241      !!----------------------------------------------------------------------------------
242      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk
243      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m]
244      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst  ! SST                                   [K]
245      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K]
246      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg]
247      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg]
248      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub   ! bulk wind speed                     [m/s]
249      !
250      INTEGER  ::   ji, jj                                ! dummy loop indices
251      REAL(wp) ::   zqa, zta, zgamma, zdth_v, ztv, zsstv  ! local scalars
252      !!-------------------------------------------------------------------
253      !
254      DO jj = 1, jpj
255         DO ji = 1, jpi
256            !
257            zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer...
258            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
259            zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer
260            zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer
261            !
262            zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!)
263            !
264            zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature"
265            !
266            ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer
267            !
268            Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk
269            !
270         END DO
271      END DO
272   END FUNCTION Ri_bulk
273
274
275   FUNCTION q_sat( ptak, pslp )
276      !!----------------------------------------------------------------------------------
277      !!                           ***  FUNCTION q_sat  ***
278      !!
279      !! ** Purpose : Specific humidity at saturation in [kg/kg]
280      !!              Based on accurate estimate of "e_sat"
281      !!              aka saturation water vapor (Goff, 1957)
282      !!
283      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
284      !!----------------------------------------------------------------------------------
285      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
286      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
287      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
288      !
289      INTEGER  ::   ji, jj         ! dummy loop indices
290      REAL(wp) ::   zta, ze_sat, ztmp   ! local scalar
291      !!----------------------------------------------------------------------------------
292      !
293      DO jj = 1, jpj
294         DO ji = 1, jpi
295            !
296            zta = MAX( ptak(ji,jj) , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
297            ztmp = rt0 / zta
298            !
299            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
300            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        &
301               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  &
302               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
303            !
304            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
305            !
306         END DO
307      END DO
308      !
309   END FUNCTION q_sat
310
311
312   !!======================================================================
313END MODULE sbcblk_phy
Note: See TracBrowser for help on using the repository browser.