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

Last change on this file since 4714 was 4714, checked in by evignon, 8 months ago

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