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

Last change on this file since 5006 was 4884, checked in by evignon, 9 months ago

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