source: LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 @ 4663

Last change on this file since 4663 was 4663, checked in by evignon, 2 years ago

ajout d'une longueur de melange sensible au cisaillement de vent dans atke
File size: 16.5 KB
RevLine 
[4449]1module atke_exchange_coeff_mod
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
[4653]28USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual
29USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd, rv, atke_ok_vdiff
[4663]30USE atke_turbulence_ini_mod, 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)
184             !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010
185             l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) &
186                                   *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0)
[4631]187          ENDIF
188      ENDDO
189   ENDDO
190
[4632]191ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
[4663]192! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
[4632]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)
[4663]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             !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010   
201             l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) &
202                                   *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0)
203
[4632]204          ENDIF
205      ENDDO
206   ENDDO
207
208
209
[4631]210ELSE
[4632]211! default: neglect effect of local stratification and shear
[4631]212
213   DO ilay=2,nlay+1
214      DO igrid=1,ngrid
215          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
216      ENDDO
217
218   ENDDO
219ENDIF
220
221
[4644]222! Computing the TKE k>=2:
223!========================
[4449]224IF (iflag_atke == 0) THEN
225
[4644]226! stationary solution (dtke/dt=0)
227
[4631]228   DO ilay=2,nlay
[4449]229        DO igrid=1,ngrid
230            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
[4631]231            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4449]232        ENDDO
233    ENDDO
234
[4631]235ELSE IF (iflag_atke == 1) THEN
236
237! full implicit scheme resolved with a second order polynomial equation
238
239    DO ilay=2,nlay
240        DO igrid=1,ngrid
[4644]241           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
[4631]242           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
243                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) &
244                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
245           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
246           tke(igrid,ilay)=0.5*(qq**2)
247        ENDDO
248    ENDDO
249
[4644]250
[4631]251ELSE IF (iflag_atke == 2) THEN
252
[4644]253! semi implicit scheme when l does not depend on tke
254! positive-guaranteed if pr slope in stable condition >1
255
256   DO ilay=2,nlay
257        DO igrid=1,ngrid
258           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
259           qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      &
260               *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
261               /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
262           tke(igrid,ilay)=0.5*(qq**2)
263        ENDDO
264    ENDDO
265
266
267ELSE IF (iflag_atke == 3) THEN
268! numerical resolution adapted from that in MAR (Deleersnijder 1992)
269! positively defined by construction
270
[4631]271    DO ilay=2,nlay
272        DO igrid=1,ngrid
[4644]273           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
274           IF (Ri(igrid,ilay) .LT. 0.) THEN
275              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
276              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
[4631]277           ELSE
[4644]278              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
279                      l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
280              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
[4631]281           ENDIF
[4644]282           qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
[4631]283           tke(igrid,ilay)=0.5*(qq**2)
284        ENDDO
285    ENDDO
286
[4644]287ELSE IF (iflag_atke == 4) THEN
288! semi implicit scheme from Arpege (V. Masson methodology with
289! Taylor expansion of the dissipation term)
[4631]290    DO ilay=2,nlay
291        DO igrid=1,ngrid
[4644]292           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
293           qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
294             +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
295             /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
296           qq=max(0.,qq)
297           tke(igrid,ilay)=0.5*(qq**2)
[4631]298        ENDDO
299    ENDDO
300
301
302ELSE
[4463]303   call abort_physic("atke_compute_km_kh", &
[4631]304        'numerical treatment of TKE not possible yet', 1)
[4449]305
306END IF
307
[4644]308! We impose a 0 tke at nlay+1
309!==============================
[4449]310
[4644]311DO igrid=1,ngrid
312 tke(igrid,nlay+1)=0.
313END DO
314
315
316! Calculation of surface TKE (k=1)
317!=================================
318! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
319DO igrid=1,ngrid
320 ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
321 tke(igrid,1)=ctkes*(ustar**2)
322END DO
323
324
325! vertical diffusion of TKE
326!==========================
[4653]327IF (atke_ok_vdiff) THEN
[4644]328   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
329ENDIF
330
331
[4449]332! Computing eddy diffusivity coefficients:
333!========================================
[4478]334DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
[4449]335    DO igrid=1,ngrid
[4478]336        ! we add the molecular viscosity to Km,h
337        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
338        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
[4449]339    END DO
340END DO
341
[4478]342! for output:
343!===========
344Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
345Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
[4449]346
347end subroutine atke_compute_km_kh
348
[4644]349!===============================================================================================
350subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
[4449]351
[4644]352! routine that computes the vertical diffusion of TKE by the turbulence
[4653]353! using an implicit resolution (See note by Dufresne and Ghattas (2009))
[4644]354! E Vignon, July 2023
355
356USE atke_turbulence_ini_mod, ONLY : rd, cke, viscom
357
358
359INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
360INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
361
362REAL, INTENT(IN)    :: dtime ! physics time step (s)
363REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
364REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
365REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
366REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
367REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: l_exchange     ! mixing length at interfaces between layers
368REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
369
370REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
371
372
373
374INTEGER                                       :: igrid,ilay
375REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
376REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
377REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
378REAL                                          :: gammak,Kem,KKb,KKt
379
380
381! Few initialisations
382CCK(:,:)=0.
383DDK(:,:)=0.
384dtke(:,:)=0.
385
386
387! Eddy diffusivity for TKE
388
389DO ilay=2,nlay
390    DO igrid=1,ngrid
391       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke
392    ENDDO
393ENDDO
394! at the top of the atmosphere set to 0
395Ke(:,nlay+1)=0.
396! at the surface, set it equal to that at the first model level
397Ke(:,1)=Ke(:,2)
398
399
400! calculate intermediary variables
401
402DO ilay=2,nlay
403    DO igrid=1,ngrid
404    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))   
405    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
406    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
407    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
408    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
409    ak(igrid,ilay)=-gammak*dtime*KKb
410    ck(igrid,ilay)=-gammak*dtime*KKt
411    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
412    ENDDO
413ENDDO
414
415! calculate CCK and DDK coefficients
416! downhill phase
417
418DO igrid=1,ngrid
419  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
420  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
421ENDDO
422
423
424DO ilay=nlay-1,2,-1
425    DO igrid=1,ngrid
426        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
427                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
428        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
429    ENDDO
430ENDDO
431
432! calculate TKE
433! uphill phase
434
435DO ilay=2,nlay+1
436    DO igrid=1,ngrid
437        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
438    ENDDO
439ENDDO
440
441! update TKE
442tke(:,:)=tke(:,:)+dtke(:,:)
443
444
445end subroutine atke_vdiff_tke
446
447
448
[4449]449end module atke_exchange_coeff_mod
Note: See TracBrowser for help on using the repository browser.