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_algo_ncar.F90 in NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 13284

Last change on this file since 13284 was 13284, checked in by smasson, 4 years ago

4.0-HEAD: merge 4.0-HEAD_r12713_clem_dan_fixcpl into 4.0-HEAD

  • Property svn:keywords set to Id
File size: 21.0 KB
Line 
1MODULE sbcblk_algo_ncar
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_ncar  ***
4   !! Computes:
5   !!   * bulk transfer coefficients C_D, C_E and C_H
6   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
7   !!   * the effective bulk wind speed at 10m U_blk
8   !!   => all these are used in bulk formulas in sbcblk.F90
9   !!
10   !!    Using the bulk formulation/param. of Large & Yeager 2008
11   !!
12   !!       Routine turb_ncar maintained and developed in AeroBulk
13   !!                     (https://github.com/brodeau/aerobulk/)
14   !!
15   !!                         L. Brodeau, 2020
16   !!=====================================================================
17   !! History :  4.0  !  2020-06  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   turb_ncar  : computes the bulk turbulent transfer coefficients
22   !!                   adjusts t_air and q_air from zt to zu m
23   !!                   returns the effective bulk wind speed at 10m
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE phycst          ! physical constants
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE sbcwave, ONLY   :  cdn_wave ! wave module
30#if defined key_si3 || defined key_cice
31   USE sbc_ice         ! Surface boundary condition: ice fields
32#endif
33   !
34   USE iom             ! I/O manager library
35   USE lib_mpp         ! distribued memory computing library
36   USE in_out_manager  ! I/O manager
37   USE prtctl          ! Print control
38   USE lib_fortran     ! to use key_nosignedzero
39
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC :: TURB_NCAR   ! called by sbcblk.F90
45
46   !                              ! NCAR own values for given constants:
47   REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature...
48
49   INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
50
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
55      &                  Cd, Ch, Ce, t_zu, q_zu, Ub,         &
56      &                  CdN, ChN, CeN                       )
57      !!----------------------------------------------------------------------
58      !!                      ***  ROUTINE  turb_ncar  ***
59      !!
60      !! ** Purpose :   Computes turbulent transfert coefficients of surface
61      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
62      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
63      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
64      !!
65      !! INPUT :
66      !! -------
67      !!    *  zt   : height for temperature and spec. hum. of air            [m]
68      !!    *  zu   : height for wind speed (usually 10m)                     [m]
69      !!    *  sst  : bulk SST                                                [K]
70      !!    *  t_zt : potential air temperature at zt                         [K]
71      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
72      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
73      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
74      !!
75      !! OUTPUT :
76      !! --------
77      !!    *  Cd     : drag coefficient
78      !!    *  Ch     : sensible heat coefficient
79      !!    *  Ce     : evaporation coefficient
80      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
81      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
82      !!    *  Ub  : bulk wind speed at zu                                 [m/s]
83      !!
84      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
85      !!----------------------------------------------------------------------------------
86      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
87      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
88      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
89      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
90      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
92      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
93      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
94      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
95      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
96      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
97      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ub    ! bulk wind speed at zu                     [m/s]
99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   CdN, ChN, CeN ! neutral transfer coefficients
100      !
101      INTEGER :: j_itt
102      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
103      !
104      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient
105      REAL(wp), DIMENSION(jpi,jpj) ::   sqrtCdN10   ! root square of Cd_n10
106      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
107      REAL(wp), DIMENSION(jpi,jpj) ::   zpsi_h_u
108      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
109      REAL(wp), DIMENSION(jpi,jpj) ::   sqrtCd
110      !!----------------------------------------------------------------------------------
111
112      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp )
113
114      Ub = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
115
116      !! Neutral drag coefficient at zu:
117      IF( ln_cdgw ) THEN      ! wave drag case
118         CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp )
119      ELSE
120         CdN = CD_N10_NCAR( Ub )
121      ENDIF
122      sqrtCdN10 = SQRT( CdN )
123
124      !! Initializing transf. coeff. with their first guess neutral equivalents :
125      Cd = CdN
126      Ce = CE_N10_NCAR( sqrtCdN10 )
127      ztmp0 = 0.5_wp + SIGN(0.5_wp, virt_temp(t_zt, q_zt) - virt_temp(sst, ssq)) ! we guess stability based on delta of virt. pot. temp.
128      Ch = CH_N10_NCAR( sqrtCdN10 , ztmp0 )
129      sqrtCd = sqrtCdN10
130
131      !! Initializing values at z_u with z_t values:
132      t_zu = t_zt
133      q_zu = q_zt
134
135      !! ITERATION BLOCK
136      DO j_itt = 1, nb_itt
137         !
138         ztmp1 = t_zu - sst   ! Updating air/sea differences
139         ztmp2 = q_zu - ssq
140
141         ! Updating turbulent scales :   (L&Y 2004 Eq. (7))
142         ztmp0 = sqrtCd*Ub          ! u*
143         ztmp1 = Ch/sqrtCd*ztmp1    ! theta*
144         ztmp2 = Ce/sqrtCd*ztmp2    ! q*
145
146         ! Estimate the inverse of Obukov length (1/L) at height zu:
147         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
148
149         !! Stability parameters :
150         zeta_u   = zu*ztmp0
151         zeta_u   = sign( min(abs(zeta_u),10._wp), zeta_u )
152
153         !! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b-9c))
154         IF( .NOT. l_zt_equal_zu ) THEN
155            ztmp0 = zt*ztmp0 ! zeta_t !
156            ztmp0 = SIGN( MIN(ABS(ztmp0),10._wp), ztmp0 )  ! Temporaty array ztmp0 == zeta_t !!!
157            ztmp0 = LOG(zt/zu) + psi_h_ncar(zeta_u) - psi_h_ncar(ztmp0)                   ! ztmp0 just used as temp array again!
158            t_zu = t_zt - ztmp1/vkarmn*ztmp0    ! ztmp1 is still theta*  L&Y 2004 Eq. (9b)
159            !!
160            q_zu = q_zt - ztmp2/vkarmn*ztmp0    ! ztmp2 is still q*      L&Y 2004 Eq. (9c)
161            q_zu = MAX(0._wp, q_zu)
162         END IF
163
164         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)...
165         !   In very rare low-wind conditions, the old way of estimating the
166         !   neutral wind speed at 10m leads to a negative value that causes the code
167         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
168         ztmp2 = psi_m_ncar(zeta_u)
169         ztmp0 = MAX( 0.25_wp , UN10_from_CD(zu, Ub, Cd, ppsi=ztmp2) ) ! U_n10 (ztmp2 == psi_m_ncar(zeta_u))
170
171         IF( ln_cdgw ) THEN      ! wave drag case
172            CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp )
173         ELSE
174            CdN   = CD_N10_NCAR(ztmp0)                                       ! Cd_n10
175         END IF
176         sqrtCdN10 = SQRT(CdN)
177
178         !! Update of transfer coefficients:
179         ztmp1  = 1._wp + sqrtCdN10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 Eq. (10a) (ztmp2 == psi_m(zeta_u))
180         Cd     = MAX( CdN / ( ztmp1*ztmp1 ) , 0.1E-3_wp )
181         sqrtCd = SQRT( Cd )
182
183         ztmp0  = ( LOG(zu/10._wp) - psi_h_ncar(zeta_u) ) / vkarmn / sqrtCdN10
184         ztmp2  = sqrtCd / sqrtCdN10
185
186         ztmp1  = 0.5_wp + sign(0.5_wp,zeta_u) ! stability flag
187         ChN    = CH_N10_NCAR( sqrtCdN10 , ztmp1 )
188         ztmp1  = 1._wp + ChN*ztmp0
189         Ch     = MAX( ChN*ztmp2 / ztmp1 , 0.1E-3_wp )   ! L&Y 2004 Eq. (10b)
190
191         CeN = CE_N10_NCAR( sqrtCdN10 )
192         ztmp1  = 1._wp + CeN*ztmp0
193         Ce     = MAX( CeN*ztmp2 / ztmp1 , 0.1E-3_wp )  ! L&Y 2004 Eq. (10c)
194
195      END DO !DO j_itt = 1, nb_itt
196
197   END SUBROUTINE turb_ncar
198
199
200   FUNCTION CD_N10_NCAR( pw10 )
201      !!----------------------------------------------------------------------------------
202      !! Estimate of the neutral drag coefficient at 10m as a function
203      !! of neutral wind  speed at 10m
204      !!
205      !! Origin: Large & Yeager 2008, Eq. (11)
206      !!
207      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
208      !!----------------------------------------------------------------------------------
209      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
210      REAL(wp), DIMENSION(jpi,jpj)             :: CD_N10_NCAR
211      !
212      INTEGER  ::     ji, jj     ! dummy loop indices
213      REAL(wp) :: zgt33, zw, zw6 ! local scalars
214      !!----------------------------------------------------------------------------------
215      !
216      DO jj = 1, jpj
217         DO ji = 1, jpi
218            !
219            zw  = pw10(ji,jj)
220            zw6 = zw*zw*zw
221            zw6 = zw6*zw6
222            !
223            ! When wind speed > 33 m/s => Cyclone conditions => special treatment
224            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
225            !
226            CD_N10_NCAR(ji,jj) = 1.e-3_wp * ( &
227               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s
228               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s
229            !
230            CD_N10_NCAR(ji,jj) = MAX( CD_N10_NCAR(ji,jj), 0.1E-3_wp )
231            !
232         END DO
233      END DO
234      !
235   END FUNCTION CD_N10_NCAR
236
237
238
239   FUNCTION CH_N10_NCAR( psqrtcdn10 , pstab )
240      !!----------------------------------------------------------------------------------
241      !! Estimate of the neutral heat transfer coefficient at 10m      !!
242      !! Origin: Large & Yeager 2008, Eq. (9) and (12)
243
244      !!----------------------------------------------------------------------------------
245      REAL(wp), DIMENSION(jpi,jpj)             :: ch_n10_ncar
246      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
247      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab      ! stable ABL => 1 / unstable ABL => 0
248      !!----------------------------------------------------------------------------------
249      !
250      ch_n10_ncar = MAX( 1.e-3_wp * psqrtcdn10*( 18._wp*pstab + 32.7_wp*(1._wp - pstab) )  , 0.1E-3_wp )   ! Eq. (9) & (12) Large & Yeager, 2008
251      !
252   END FUNCTION CH_N10_NCAR
253
254   FUNCTION CE_N10_NCAR( psqrtcdn10 )
255      !!----------------------------------------------------------------------------------
256      !! Estimate of the neutral heat transfer coefficient at 10m      !!
257      !! Origin: Large & Yeager 2008, Eq. (9) and (13)
258      !!----------------------------------------------------------------------------------
259      REAL(wp), DIMENSION(jpi,jpj)             :: ce_n10_ncar
260      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
261      !!----------------------------------------------------------------------------------
262      ce_n10_ncar = MAX( 1.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , 0.1E-3_wp )
263      !
264   END FUNCTION CE_N10_NCAR
265
266
267   FUNCTION psi_m_ncar( pzeta )
268      !!----------------------------------------------------------------------------------
269      !! Universal profile stability function for momentum
270      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
271      !!
272      !! pzeta : stability paramenter, z/L where z is altitude measurement
273      !!         and L is M-O length
274      !!
275      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
276      !!----------------------------------------------------------------------------------
277      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar
278      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
279      !
280      INTEGER  ::   ji, jj    ! dummy loop indices
281      REAL(wp) :: zzeta, zx2, zx, zpsi_unst, zpsi_stab,  zstab   ! local scalars
282      !!----------------------------------------------------------------------------------
283      DO jj = 1, jpj
284         DO ji = 1, jpi
285
286            zzeta = pzeta(ji,jj)
287            !
288            zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 - 16z)^0.5
289            zx2 = MAX( zx2 , 1._wp )
290            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25
291            zpsi_unst = 2._wp*LOG( (1._wp + zx )*0.5_wp )   &
292               &            + LOG( (1._wp + zx2)*0.5_wp )   &
293               &          - 2._wp*ATAN(zx) + rpi*0.5_wp
294            !
295            zpsi_stab = -5._wp*zzeta
296            !
297            zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1
298            !
299            psi_m_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable
300               &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable
301            !
302         END DO
303      END DO
304   END FUNCTION psi_m_ncar
305
306
307   FUNCTION psi_h_ncar( pzeta )
308      !!----------------------------------------------------------------------------------
309      !! Universal profile stability function for temperature and humidity
310      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
311      !!
312      !! pzeta : stability paramenter, z/L where z is altitude measurement
313      !!         and L is M-O length
314      !!
315      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
316      !!----------------------------------------------------------------------------------
317      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar
318      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
319      !
320      INTEGER  ::   ji, jj     ! dummy loop indices
321      REAL(wp) :: zzeta, zx2, zpsi_unst, zpsi_stab, zstab  ! local scalars
322      !!----------------------------------------------------------------------------------
323      !
324      DO jj = 1, jpj
325         DO ji = 1, jpi
326            !
327            zzeta = pzeta(ji,jj)
328            !
329            zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 -16z)^0.5
330            zx2 = MAX( zx2 , 1._wp )
331            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
332            !
333            zpsi_stab = -5._wp*zzeta
334            !
335            zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1
336            !
337            psi_h_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable
338               &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable
339            !
340         END DO
341      END DO
342   END FUNCTION psi_h_ncar
343
344
345
346
347   FUNCTION UN10_from_CD( pzu, pUb, pCd, ppsi )
348      !!----------------------------------------------------------------------------------
349      !!  Provides the neutral-stability wind speed at 10 m
350      !!----------------------------------------------------------------------------------
351      REAL(wp), DIMENSION(jpi,jpj)             :: UN10_from_CD  !: [m/s]
352      REAL(wp),                     INTENT(in) :: pzu  !: measurement heigh of bulk wind speed
353      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb  !: bulk wind speed at height pzu m   [m/s]
354      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd  !: drag coefficient
355      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
356      !!----------------------------------------------------------------------------------
357      !! Reminder: UN10 = u*/vkarmn * log(10/z0)
358      !!     and: u* = sqrt(Cd) * Ub
359      !!                                  u*/vkarmn * log(   10   /       z0    )
360      UN10_from_CD(:,:) = SQRT(pCd(:,:))*pUb/vkarmn * LOG( 10._wp / z0_from_Cd( pzu, pCd(:,:), ppsi=ppsi(:,:) ) )
361      !!
362   END FUNCTION UN10_from_CD
363
364
365   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
366      !!------------------------------------------------------------------------
367      !!
368      !! Evaluates the 1./(Obukhov length) from air temperature,
369      !! air specific humidity, and frictional scales u*, t* and q*
370      !!
371      !! Author: L. Brodeau, June 2019 / AeroBulk
372      !!         (https://github.com/brodeau/aerobulk/)
373      !!------------------------------------------------------------------------
374      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L     !: 1./(Obukhov length) [m^-1]
375      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha         !: reference potential temperature of air [K]
376      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa          !: reference specific humidity of air   [kg/kg]
377      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus          !: u*: friction velocity [m/s]
378      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts, pqs     !: \theta* and q* friction aka turb. scales for temp. and spec. hum.
379      !
380      INTEGER  ::   ji, jj         ! dummy loop indices
381      REAL(wp) ::     zqa          ! local scalar
382      !!-------------------------------------------------------------------
383      !
384      DO jj = 1, jpj
385         DO ji = 1, jpi
386            !
387            zqa = (1._wp + rctv0*pqa(ji,jj))
388            !
389            ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
390            !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
391            !                      or
392            !  b/  -u* [ theta*              + 0.61 theta q* ]
393            !
394            One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
395               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
396            !
397         END DO
398      END DO
399      !
400      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
401      !
402   END FUNCTION One_on_L
403
404
405   FUNCTION z0_from_Cd( pzu, pCd,  ppsi )
406      REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd        !: roughness length [m]
407      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m]
408      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd   !: (neutral or non-neutral) drag coefficient []
409      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
410      !!
411      !! If pCd is the NEUTRAL-STABILITY drag coefficient then ppsi must be 0 or not given
412      !! If pCd is the drag coefficient (in stable or unstable conditions) then pssi must be provided
413      !!----------------------------------------------------------------------------------
414      IF ( PRESENT(ppsi) ) THEN
415         !! Cd provided is the actual Cd (not the neutral-stability CdN) :
416         z0_from_Cd = pzu * EXP( - ( vkarmn/SQRT(pCd(:,:)) + ppsi(:,:) ) ) !LB: ok, double-checked!
417      ELSE
418         !! Cd provided is the neutral-stability Cd, aka CdN :
419         z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) )            !LB: ok, double-checked!
420      END IF
421   END FUNCTION z0_from_Cd
422
423   FUNCTION virt_temp( pta, pqa )
424      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp !: virtual temperature [K]
425      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K]
426      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air   [kg/kg]
427      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
428   END FUNCTION virt_temp
429
430   !!======================================================================
431END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.