source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_atke_exchange_coeff.F90 @ 5441

Last change on this file since 5441 was 5119, checked in by abarral, 5 months ago

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