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/branches/2020/temporary_r4_trunk/src/OCE/SBC – NEMO

source: NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 13470

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

r4_trunk: second change of DO loops for routines to be merged, see #2523
  • Property svn:keywords set to Id
File size: 20.7 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_2D( 1, 1, 1, 1 )
217         !
218         zw  = pw10(ji,jj)
219         zw6 = zw*zw*zw
220         zw6 = zw6*zw6
221         !
222         ! When wind speed > 33 m/s => Cyclone conditions => special treatment
223         zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
224         !
225         CD_N10_NCAR(ji,jj) = 1.e-3_wp * ( &
226            &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s
227            &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s
228         !
229         CD_N10_NCAR(ji,jj) = MAX( CD_N10_NCAR(ji,jj), 0.1E-3_wp )
230         !
231      END_2D
232      !
233   END FUNCTION CD_N10_NCAR
234
235
236
237   FUNCTION CH_N10_NCAR( psqrtcdn10 , pstab )
238      !!----------------------------------------------------------------------------------
239      !! Estimate of the neutral heat transfer coefficient at 10m      !!
240      !! Origin: Large & Yeager 2008, Eq. (9) and (12)
241
242      !!----------------------------------------------------------------------------------
243      REAL(wp), DIMENSION(jpi,jpj)             :: ch_n10_ncar
244      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
245      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab      ! stable ABL => 1 / unstable ABL => 0
246      !!----------------------------------------------------------------------------------
247      !
248      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
249      !
250   END FUNCTION CH_N10_NCAR
251
252   FUNCTION CE_N10_NCAR( psqrtcdn10 )
253      !!----------------------------------------------------------------------------------
254      !! Estimate of the neutral heat transfer coefficient at 10m      !!
255      !! Origin: Large & Yeager 2008, Eq. (9) and (13)
256      !!----------------------------------------------------------------------------------
257      REAL(wp), DIMENSION(jpi,jpj)             :: ce_n10_ncar
258      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
259      !!----------------------------------------------------------------------------------
260      ce_n10_ncar = MAX( 1.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , 0.1E-3_wp )
261      !
262   END FUNCTION CE_N10_NCAR
263
264
265   FUNCTION psi_m_ncar( pzeta )
266      !!----------------------------------------------------------------------------------
267      !! Universal profile stability function for momentum
268      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
269      !!
270      !! pzeta : stability paramenter, z/L where z is altitude measurement
271      !!         and L is M-O length
272      !!
273      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
274      !!----------------------------------------------------------------------------------
275      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar
276      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
277      !
278      INTEGER  ::   ji, jj    ! dummy loop indices
279      REAL(wp) :: zzeta, zx2, zx, zpsi_unst, zpsi_stab,  zstab   ! local scalars
280      !!----------------------------------------------------------------------------------
281      DO_2D( 1, 1, 1, 1 )
282
283         zzeta = pzeta(ji,jj)
284         !
285         zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 - 16z)^0.5
286         zx2 = MAX( zx2 , 1._wp )
287         zx  = SQRT(zx2)                          ! (1 - 16z)^0.25
288         zpsi_unst = 2._wp*LOG( (1._wp + zx )*0.5_wp )   &
289            &            + LOG( (1._wp + zx2)*0.5_wp )   &
290            &          - 2._wp*ATAN(zx) + rpi*0.5_wp
291         !
292         zpsi_stab = -5._wp*zzeta
293         !
294         zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1
295         !
296         psi_m_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable
297            &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable
298         !
299      END_2D
300   END FUNCTION psi_m_ncar
301
302
303   FUNCTION psi_h_ncar( pzeta )
304      !!----------------------------------------------------------------------------------
305      !! Universal profile stability function for temperature and humidity
306      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
307      !!
308      !! pzeta : stability paramenter, z/L where z is altitude measurement
309      !!         and L is M-O length
310      !!
311      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
312      !!----------------------------------------------------------------------------------
313      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar
314      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
315      !
316      INTEGER  ::   ji, jj     ! dummy loop indices
317      REAL(wp) :: zzeta, zx2, zpsi_unst, zpsi_stab, zstab  ! local scalars
318      !!----------------------------------------------------------------------------------
319      !
320      DO_2D( 1, 1, 1, 1 )
321         !
322         zzeta = pzeta(ji,jj)
323         !
324         zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 -16z)^0.5
325         zx2 = MAX( zx2 , 1._wp )
326         zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
327         !
328         zpsi_stab = -5._wp*zzeta
329         !
330         zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1
331         !
332         psi_h_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable
333            &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable
334         !
335      END_2D
336   END FUNCTION psi_h_ncar
337
338
339
340
341   FUNCTION UN10_from_CD( pzu, pUb, pCd, ppsi )
342      !!----------------------------------------------------------------------------------
343      !!  Provides the neutral-stability wind speed at 10 m
344      !!----------------------------------------------------------------------------------
345      REAL(wp), DIMENSION(jpi,jpj)             :: UN10_from_CD  !: [m/s]
346      REAL(wp),                     INTENT(in) :: pzu  !: measurement heigh of bulk wind speed
347      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb  !: bulk wind speed at height pzu m   [m/s]
348      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd  !: drag coefficient
349      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
350      !!----------------------------------------------------------------------------------
351      !! Reminder: UN10 = u*/vkarmn * log(10/z0)
352      !!     and: u* = sqrt(Cd) * Ub
353      !!                                  u*/vkarmn * log(   10   /       z0    )
354      UN10_from_CD(:,:) = SQRT(pCd(:,:))*pUb/vkarmn * LOG( 10._wp / z0_from_Cd( pzu, pCd(:,:), ppsi=ppsi(:,:) ) )
355      !!
356   END FUNCTION UN10_from_CD
357
358
359   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
360      !!------------------------------------------------------------------------
361      !!
362      !! Evaluates the 1./(Obukhov length) from air temperature,
363      !! air specific humidity, and frictional scales u*, t* and q*
364      !!
365      !! Author: L. Brodeau, June 2019 / AeroBulk
366      !!         (https://github.com/brodeau/aerobulk/)
367      !!------------------------------------------------------------------------
368      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L     !: 1./(Obukhov length) [m^-1]
369      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha         !: reference potential temperature of air [K]
370      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa          !: reference specific humidity of air   [kg/kg]
371      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus          !: u*: friction velocity [m/s]
372      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts, pqs     !: \theta* and q* friction aka turb. scales for temp. and spec. hum.
373      !
374      INTEGER  ::   ji, jj         ! dummy loop indices
375      REAL(wp) ::     zqa          ! local scalar
376      !!-------------------------------------------------------------------
377      !
378      DO_2D( 1, 1, 1, 1 )
379         !
380         zqa = (1._wp + rctv0*pqa(ji,jj))
381         !
382         ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
383         !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
384         !                      or
385         !  b/  -u* [ theta*              + 0.61 theta q* ]
386         !
387         One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
388            &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
389         !
390      END_2D
391      !
392      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
393      !
394   END FUNCTION One_on_L
395
396
397   FUNCTION z0_from_Cd( pzu, pCd,  ppsi )
398      REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd        !: roughness length [m]
399      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m]
400      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd   !: (neutral or non-neutral) drag coefficient []
401      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
402      !!
403      !! If pCd is the NEUTRAL-STABILITY drag coefficient then ppsi must be 0 or not given
404      !! If pCd is the drag coefficient (in stable or unstable conditions) then pssi must be provided
405      !!----------------------------------------------------------------------------------
406      IF ( PRESENT(ppsi) ) THEN
407         !! Cd provided is the actual Cd (not the neutral-stability CdN) :
408         z0_from_Cd = pzu * EXP( - ( vkarmn/SQRT(pCd(:,:)) + ppsi(:,:) ) ) !LB: ok, double-checked!
409      ELSE
410         !! Cd provided is the neutral-stability Cd, aka CdN :
411         z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) )            !LB: ok, double-checked!
412      END IF
413   END FUNCTION z0_from_Cd
414
415   FUNCTION virt_temp( pta, pqa )
416      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp !: virtual temperature [K]
417      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K]
418      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air   [kg/kg]
419      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
420   END FUNCTION virt_temp
421
422   !!======================================================================
423END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.