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

Last change on this file since 5173 was 5119, checked in by abarral, 12 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.