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_phymbl.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_phymbl.F90 @ 11111

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

LB: skin temperature now used by ECMWF bulk algorithm, and can be xios-ed ; new module "SBC/sbcblk_skin.F90" ; all physical constants defined in SBC now moved to "DOM/phycst.F90"; all air-properties and other ABL functions gathered in a new module: "SBC/sbcblk_phymbl.F90".
File size: 9.0 KB
Line 
1MODULE sbcblk_phymbl
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_phymbl  ***
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   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
11   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
12   !!   q_sat         : saturation humidity as a function of SLP and temperature
13   !!   gamma_moist   :
14   !!   L_vap         : latent heat of vaporization of water as a function of temperature
15   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature
16   
17   USE oce            ! ocean dynamics and tracers
18   USE dom_oce        ! ocean space and time domain
19   USE phycst         ! physical constants
20   !USE sbc_oce        ! Surface boundary condition: ocean fields
21   !USE sbcdcy         ! surface boundary condition: diurnal cycle
22   !USE sbcwave , ONLY :   cdn_wave ! wave module
23   !USE sbc_ice        ! Surface boundary condition: ice fields
24   USE lib_fortran    ! to use key_nosignedzero
25   !
26   !USE iom            ! I/O manager library
27   !USE in_out_manager ! I/O manager
28   USE lib_mpp        ! distribued memory computing library
29   !USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
30   !USE prtctl         ! Print control
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC rho_air
36   PUBLIC cp_air
37   PUBLIC q_sat
38   PUBLIC gamma_moist
39   PUBLIC L_vap
40   PUBLIC visc_air
41
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   FUNCTION rho_air( ptak, pqa, pslp )
50      !!-------------------------------------------------------------------------------
51      !!                           ***  FUNCTION rho_air  ***
52      !!
53      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
54      !!
55      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
56      !!-------------------------------------------------------------------------------
57      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
58      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
59      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
60      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
61      !!-------------------------------------------------------------------------------
62      !
63      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
64      !
65   END FUNCTION rho_air
66
67
68   FUNCTION cp_air( pqa )
69      !!-------------------------------------------------------------------------------
70      !!                           ***  FUNCTION cp_air  ***
71      !!
72      !! ** Purpose : Compute specific heat (Cp) of moist air
73      !!
74      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
75      !!-------------------------------------------------------------------------------
76      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
77      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
78      !!-------------------------------------------------------------------------------
79      !
80      cp_air = rCp_dry + rCp_vap * pqa
81      !
82   END FUNCTION cp_air
83
84
85   FUNCTION q_sat( ptak, pslp )
86      !!----------------------------------------------------------------------------------
87      !!                           ***  FUNCTION q_sat  ***
88      !!
89      !! ** Purpose : Specific humidity at saturation in [kg/kg]
90      !!              Based on accurate estimate of "e_sat"
91      !!              aka saturation water vapor (Goff, 1957)
92      !!
93      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
94      !!----------------------------------------------------------------------------------
95      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
96      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
97      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
98      !
99      INTEGER  ::   ji, jj         ! dummy loop indices
100      REAL(wp) ::   ze_sat, ztmp   ! local scalar
101      !!----------------------------------------------------------------------------------
102      !
103      DO jj = 1, jpj
104         DO ji = 1, jpi
105            !
106            ztmp = rt0 / ptak(ji,jj)
107            !
108            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
109            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
110               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
111               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
112            !
113            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]
114            !
115         END DO
116      END DO
117      !
118   END FUNCTION q_sat
119
120
121   FUNCTION gamma_moist( ptak, pqa )
122      !!----------------------------------------------------------------------------------
123      !!                           ***  FUNCTION gamma_moist  ***
124      !!
125      !! ** Purpose : Compute the moist adiabatic lapse-rate.
126      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
127      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
128      !!
129      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
130      !!----------------------------------------------------------------------------------
131      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
132      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
133      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
134      !
135      INTEGER  ::   ji, jj         ! dummy loop indices
136      REAL(wp) :: zrv, ziRT        ! local scalar
137      !!----------------------------------------------------------------------------------
138      !
139      DO jj = 1, jpj
140         DO ji = 1, jpi
141            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
142            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
143            gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( rCp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )
144         END DO
145      END DO
146      !
147   END FUNCTION gamma_moist
148
149
150   FUNCTION L_vap( psst )
151      !!---------------------------------------------------------------------------------
152      !!                           ***  FUNCTION L_vap  ***
153      !!
154      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
155      !!
156      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
157      !!----------------------------------------------------------------------------------
158      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
159      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
160      !!----------------------------------------------------------------------------------
161      !
162      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
163      !
164   END FUNCTION L_vap
165
166
167
168   FUNCTION visc_air(ptak)
169      !!----------------------------------------------------------------------------------
170      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
171      !!
172      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
173      !!----------------------------------------------------------------------------------
174      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
175      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
176      !
177      INTEGER  ::   ji, jj      ! dummy loop indices
178      REAL(wp) ::   ztc, ztc2   ! local scalar
179      !!----------------------------------------------------------------------------------
180      !
181      DO jj = 1, jpj
182         DO ji = 1, jpi
183            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
184            ztc2 = ztc*ztc
185            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
186         END DO
187      END DO
188      !
189   END FUNCTION visc_air
190   
191
192   !!======================================================================
193END MODULE sbcblk_phymbl
Note: See TracBrowser for help on using the repository browser.