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

Last change on this file since 4799 was 4780, checked in by evignon, 13 months ago

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