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

Last change on this file since 4744 was 4714, checked in by evignon, 10 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.