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

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

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