source: LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90 @ 5087

Last change on this file since 5087 was 5039, checked in by evignon, 12 months ago

ajout des tendances de tke de la routine ATKE dans les sorties
File size: 22.7 KB
Line 
1module lmdz_atke_exchange_coeff
2
3implicit none
4
5contains
6
7subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
8                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
9                        tke,eps,tke_shear,tke_buoy,tke_trans,Km_out,Kh_out)
10
11!========================================================================
12! Routine that computes turbulent Km / Kh coefficients with a
13! 1.5 order closure scheme (TKE) with or without stationarity assumption
14!
15! This parameterization has been constructed in the framework of a
16! collective and collaborative workshop,
17! the so-called 'Atelier TKE (ATKE)' with
18! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos,
19! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon
20!
21! Main assumptions of the model :
22! (1) horizontal homogeneity (Dx=Dy=0.)
23!=======================================================================
24
25USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual, ri0, ri1
26USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes, rg, rd, rv, atke_ok_vdiff
27USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin, cn
28
29!!-------------------------------------------------------------------------------------------------------------
30! integer :: iflag_atke        ! flag that controls options in atke_compute_km_kh
31! integer :: iflag_atke_lmix   ! flag that controls the calculation of mixing length in atke
32! integer :: iflag_num_atke    ! flag that controls the numerical treatment of diffusion coeffiient calculation
33
34! logical :: atke_ok_vdiff    ! activate vertical diffusion of TKE or not
35! logical :: atke_ok_virtual  ! account for vapor for flottability
36
37! real :: kappa = 0.4         ! Von Karman constant
38! real :: cn                  ! Sm value at Ri=0
39! real :: cinf                ! Sm value at Ri=-Inf
40! real :: ri0                 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0
41! real :: ri1                 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0
42! real :: lmin                ! minimum mixing length corresponding to the Kolmogov dissipation scale  in planetary atmospheres (Chen et al 2016, JGR Atmos)
43! real :: ctkes               ! coefficient for surface TKE
44! real :: clmixshear          ! coefficient for mixing length depending on local wind shear
45
46! Tunable parameters for the ATKE scheme and their range of values
47!!-------------------------------------------------------------------------------------------------------------
48! real :: cepsilon            ! controls the value of the dissipation length scale, range [1.2 - 10]
49! real :: cke                 ! controls the value of the diffusion coefficient of TKE, range [1 - 5]
50! real :: l0                  ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75]
51! real :: clmix               ! controls the value of the mixing length in stratified conditions, range [0.1 - 2]
52! real :: ric                 ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25]
53! real :: smmin               ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1]
54! real :: pr_neut             ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1]
55! real :: pr_slope            ! linear slope of Pr with Ri in the very stable regime, range [3 - 5]
56! real :: cinffac             ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0]
57! real :: pr_asym             ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5]
58!!-------------------------------------------------------------------------------------------------------------
59
60implicit none
61
62
63! Declarations:
64!=============
65
66INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
67INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
68
69REAL, INTENT(IN)    :: dtime ! physics time step (s)
70REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
71REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
72REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp     ! temperature (K)
73REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap     ! specific humidity (kg/kg)
74REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play     ! pressure (Pa)
75REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf  ! pressure at interfaces(Pa)
76REAL, DIMENSION(ngrid), INTENT(IN)            :: cdrag_uv ! surface drag coefficient for momentum
77
78REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke      ! turbulent kinetic energy at interface between layers (m2/s2)
79
80REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: eps      ! output: TKE dissipation rate at interface between layers (m2/s3)
81REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: tke_shear! output: TKE shear production rate (m2/s3)
82REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: tke_buoy ! output: TKE buoyancy production rate (m2/s3)
83REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: tke_trans! output: TKE transport (diffusion) term (m2/s3)
84REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers (m2/s)
85REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers (m2/s)
86
87! Local variables
88REAL, DIMENSION(ngrid,nlay+1) :: Km          ! Exchange coefficient for momentum at interface between layers
89REAL, DIMENSION(ngrid,nlay+1) :: Kh          ! Exchange coefficient for heat flux at interface between layers
90REAL, DIMENSION(ngrid,nlay)   :: theta       ! Potential temperature
91REAL, DIMENSION(ngrid,nlay+1) :: l_exchange  ! Length of exchange (at interface)
92REAL, DIMENSION(ngrid,nlay+1) :: z_interf    ! Altitude at the interface
93REAL, DIMENSION(ngrid,nlay)   :: z_lay       ! Altitude of layers
94REAL, DIMENSION(ngrid,nlay)   :: dz_interf   ! distance between two consecutive interfaces
95REAL, DIMENSION(ngrid,nlay)   :: dz_lay      ! distance between two layer middles (NB: first and last are half layers)
96REAL, DIMENSION(ngrid,nlay+1) :: N2          ! square of Brunt Vaisala pulsation (at interface)
97REAL, DIMENSION(ngrid,nlay+1) :: shear2      ! square of wind shear (at interface)
98REAL, DIMENSION(ngrid,nlay+1) :: Ri          ! Richardson's number (at interface)
99REAL, DIMENSION(ngrid,nlay+1) :: Prandtl     ! Turbulent Prandtl's number (at interface)
100REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
101REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
102
103INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
104REAL    :: preff      ! reference pressure for potential temperature calculations
105REAL    :: thetam     ! mean potential temperature at interface
106REAL    :: delta      ! discriminant of the second order polynomial
107REAL    :: qq         ! tke=qq**2/2
108REAL    :: shear      ! wind shear
109REAL    :: lstrat     ! mixing length depending on local stratification
110REAL    :: taustrat   ! caracteristic timescale for turbulence in very stable conditions
111REAL    :: netloss    ! net loss term of tke
112REAL    :: netsource  ! net source term of tke
113REAL    :: ustar      ! friction velocity estimation
114REAL    :: invtau     
115REAL    :: rvap
116
117! Initializations:
118!================
119
120DO igrid=1,ngrid
121    dz_interf(igrid,1) = 0.0
122    z_interf(igrid,1) = 0.0
123END DO
124
125! Calculation of potential temperature: (if vapor -> virtual potential temperature)
126!=====================================
127
128preff=100000.
129! results should not depend on the choice of preff
130DO ilay=1,nlay
131    DO igrid = 1, ngrid
132        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
133    END DO
134END DO
135
136! account for water vapor mass for buoyancy calculation
137IF (atke_ok_virtual) THEN
138    DO ilay=1,nlay
139        DO igrid = 1, ngrid
140            rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
141            theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
142        END DO
143    END DO
144ENDIF
145
146
147! Calculation of altitude of layers' middle and bottom interfaces:
148!=================================================================
149
150DO ilay=2,nlay+1
151    DO igrid=1,ngrid
152        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
153        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
154    ENDDO
155ENDDO
156
157DO ilay=1,nlay
158    DO igrid=1,ngrid
159        z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
160    ENDDO
161ENDDO
162
163
164! Computes the gradient Richardson's number and stability functions:
165!===================================================================
166
167DO ilay=2,nlay
168    DO igrid=1,ngrid
169        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
170        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
171        N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
172        shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
173            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 )
174        Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10)
175       
176        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
177            Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
178            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut
179        ELSE ! stable cases
180            Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric))
181            ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
182            Prandtl(igrid,ilay) = pr_neut*exp(-pr_slope/pr_neut*Ri(igrid,ilay)+Ri(igrid,ilay)/pr_neut) &
183                                + Ri(igrid,ilay) * pr_slope
184            IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
185                call abort_physic("atke_compute_km_kh", &
186                'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
187            ENDIF
188        END IF
189       
190        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
191
192    ENDDO
193ENDDO
194
195
196! Computing the mixing length:
197!==============================================================
198
199
200IF (iflag_atke_lmix .EQ. 1 ) THEN
201    ! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions
202        DO ilay=2,nlay
203            DO igrid=1,ngrid
204                l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
205                IF (N2(igrid,ilay) .GT. 0.) THEN
206                    lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
207                    lstrat=max(lstrat,lmin)
208                    !Inverse interpolation, Van de Wiel et al. 2010
209                    l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
210                ENDIF
211            ENDDO
212        ENDDO
213
214    ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
215    ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
216    ! implies 2 tuning coefficients clmix and clmixshear
217    DO ilay=2,nlay
218        DO igrid=1,ngrid
219            l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
220            IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
221                lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &
222                    clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay)))
223                lstrat=max(lstrat,lmin)
224                !Inverse interpolation, Van de Wiel et al. 2010   
225                l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
226            ENDIF
227        ENDDO
228    ENDDO
229
230    ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN
231    ! add effect of wind shear on lstrat following grisogono 2010, qjrms
232    ! keeping a single tuning coefficient clmix
233    DO ilay=2,nlay
234        DO igrid=1,ngrid
235            l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
236            IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
237                lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay)))
238                lstrat=max(lstrat,lmin)
239                !Inverse interpolation, Van de Wiel et al. 2010   
240                l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
241            ENDIF
242        ENDDO
243    ENDDO
244
245ELSE
246    ! default Blackadar formulation: neglect effect of local stratification and shear
247    DO ilay=2,nlay+1
248        DO igrid=1,ngrid
249            l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
250        ENDDO
251    ENDDO
252ENDIF
253
254
255! Computing the TKE k>=2:
256!========================
257IF (iflag_atke == 0) THEN
258
259    ! stationary solution (dtke/dt=0)
260
261    DO ilay=2,nlay
262            DO igrid=1,ngrid
263                tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
264                shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
265                eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
266                tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay)
267                tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
268                                    *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
269            ENDDO
270        ENDDO
271
272ELSE IF (iflag_atke == 1) THEN
273
274    ! full implicit scheme resolved with a second order polynomial equation
275    ! default solution which shows fair convergence properties
276        DO ilay=2,nlay
277            DO igrid=1,ngrid
278            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
279            delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. &
280                    +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + &
281                    2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) &
282                    *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
283            qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
284            qq=max(0.,qq)
285            tke(igrid,ilay)=0.5*(qq**2)
286            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
287            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay)
288            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
289                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
290            ENDDO
291        ENDDO
292
293
294ELSE IF (iflag_atke == 2) THEN
295
296    ! semi implicit scheme when l does not depend on tke
297    ! positive-guaranteed if pr slope in stable condition >1
298
299    DO ilay=2,nlay
300            DO igrid=1,ngrid
301            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
302            qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      &
303                *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
304                /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
305            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
306            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
307                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
308            tke(igrid,ilay)=0.5*(qq**2)
309            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
310            ENDDO
311        ENDDO
312
313
314ELSE IF (iflag_atke == 3) THEN
315    ! numerical resolution adapted from that in MAR (Deleersnijder 1992)
316    ! positively defined by construction
317
318        DO ilay=2,nlay
319            DO igrid=1,ngrid
320            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
321            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
322            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
323            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
324                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
325            IF (Ri(igrid,ilay) .LT. 0.) THEN
326                netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
327                netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
328            ELSE
329                netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
330                        l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
331                netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
332            ENDIF
333            qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
334            tke(igrid,ilay)=0.5*(qq**2)
335            ENDDO
336        ENDDO
337
338ELSE IF (iflag_atke == 4) THEN
339    ! semi implicit scheme from Arpege (V. Masson methodology with
340    ! Taylor expansion of the dissipation term)
341        DO ilay=2,nlay
342            DO igrid=1,ngrid
343            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
344            tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
345            tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
346                                *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
347            qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
348                +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
349                /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
350            qq=max(0.,qq)
351            tke(igrid,ilay)=0.5*(qq**2)
352            eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
353            ENDDO
354        ENDDO
355
356
357ELSE
358    call abort_physic("atke_compute_km_kh", &
359        'numerical treatment of TKE not possible yet', 1)
360
361END IF
362
363! We impose a 0 tke at nlay+1
364!==============================
365
366DO igrid=1,ngrid
367    tke(igrid,nlay+1)=0.
368    eps(igrid,nlay+1)=0.
369    tke_shear(igrid,nlay+1)=0.
370    tke_buoy(igrid,nlay+1)=0.
371END DO
372
373
374! Calculation of surface TKE (k=1)
375!=================================
376! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
377DO igrid=1,ngrid
378    ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
379    tke(igrid,1)=ctkes*(ustar**2)
380    eps(igrid,1)=0. ! arbitrary as TKE is not properly defined at the surface
381    tke_shear(igrid,1)=0.
382    tke_buoy(igrid,1)=0.
383END DO
384
385
386! vertical diffusion of TKE
387!==========================
388tke_trans(:,:)=0.
389IF (atke_ok_vdiff) THEN
390    CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
391ENDIF
392
393
394! Computing eddy diffusivity coefficients:
395!========================================
396DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
397    DO igrid=1,ngrid
398        ! we add the molecular viscosity to Km,h
399        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
400        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
401    END DO
402END DO
403
404! for output:
405!===========
406Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
407Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
408
409end subroutine atke_compute_km_kh
410
411!===============================================================================================
412subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
413
414! routine that computes the vertical diffusion of TKE by the turbulence
415! using an implicit resolution (See note by Dufresne and Ghattas (2009))
416! E Vignon, July 2023
417
418USE lmdz_atke_turbulence_ini, ONLY : rd, cke, viscom
419
420
421INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
422INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
423
424REAL, INTENT(IN)    :: dtime ! physics time step (s)
425REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
426REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
427REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
428REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
429REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: l_exchange     ! mixing length at interfaces between layers
430REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
431
432REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
433REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke_trans ! turbulent kinetic energy transport term (m2/s3)
434
435
436INTEGER                                       :: igrid,ilay
437REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
438REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
439REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
440REAL                                          :: gammak,Kem,KKb,KKt
441
442
443! Few initialisations
444CCK(:,:)=0.
445DDK(:,:)=0.
446dtke(:,:)=0.
447
448
449! Eddy diffusivity for TKE
450
451DO ilay=2,nlay
452    DO igrid=1,ngrid
453        Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke
454    ENDDO
455ENDDO
456! at the top of the atmosphere set to 0
457Ke(:,nlay+1)=0.
458! at the surface, set it equal to that at the first model level
459Ke(:,1)=Ke(:,2)
460
461
462! calculate intermediary variables
463
464DO ilay=2,nlay
465    DO igrid=1,ngrid
466    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))   
467    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
468    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
469    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
470    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
471    ak(igrid,ilay)=-gammak*dtime*KKb
472    ck(igrid,ilay)=-gammak*dtime*KKt
473    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
474    ENDDO
475ENDDO
476
477! calculate CCK and DDK coefficients
478! downhill phase
479
480DO igrid=1,ngrid
481    CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
482    DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
483ENDDO
484
485
486DO ilay=nlay-1,2,-1
487    DO igrid=1,ngrid
488        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
489                        / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
490        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
491    ENDDO
492ENDDO
493
494! calculate TKE
495! uphill phase
496
497DO ilay=2,nlay+1
498    DO igrid=1,ngrid
499        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
500    ENDDO
501ENDDO
502
503! update TKE
504tke(:,:)=tke(:,:)+dtke(:,:)
505tke_trans(:,:)=dtke(:,:)/dtime
506
507
508end subroutine atke_vdiff_tke
509
510
511
512end module lmdz_atke_exchange_coeff
Note: See TracBrowser for help on using the repository browser.