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

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