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

Last change on this file since 4660 was 4653, checked in by evignon, 15 months ago

prise en compte de l'humidite pour le calcul du flux de flottabilite dans atke
+ petits ajustements
File size: 17.4 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
[4644]30USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, 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)
[4631]70LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases
[4449]71
72INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
[4481]73REAL    :: cn,Ri0,Ri1    ! parameter for Sm stability function and Prandlt
[4478]74REAL    :: preff      ! reference pressure for potential temperature calculations
75REAL    :: thetam     ! mean potential temperature at interface
[4631]76REAL    :: delta      ! discriminant of the second order polynomial
77REAL    :: qq         ! tke=qq**2/2
78REAL    :: shear      ! wind shear
79REAL    :: lstrat     ! mixing length depending on local stratification
80REAL    :: taustrat   ! caracteristic timescale for turbulence in very stable conditions
[4644]81REAL    :: netloss    ! net loss term of tke
82REAL    :: netsource  ! net source term of tke
83REAL    :: ustar      ! friction velocity estimation
[4653]84REAL    :: invtau     
85REAL    :: rvap
[4449]86
87! Initializations:
88!================
89
90DO igrid=1,ngrid
[4478]91    dz_interf(igrid,1) = 0.0
[4449]92    z_interf(igrid,1) = 0.0
93END DO
94
[4653]95! Calculation of potential temperature: (if vapor -> virtual potential temperature)
[4478]96!=====================================
[4449]97
[4478]98preff=100000.
[4653]99! results should not depend on the choice of preff
[4478]100DO ilay=1,nlay
101     DO igrid = 1, ngrid
102        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
103     END DO
104END DO
[4449]105
[4653]106! account for water vapor mass for buoyancy calculation
107IF (atke_ok_virtual) THEN
108  DO ilay=1,nlay
109     DO igrid = 1, ngrid
110        rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
111        theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
112     END DO
113  END DO
114ENDIF
[4449]115
[4478]116
117! Calculation of altitude of layers' middle and bottom interfaces:
118!=================================================================
119
[4449]120DO ilay=2,nlay+1
121    DO igrid=1,ngrid
[4478]122        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
123        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
[4449]124    ENDDO
125ENDDO
126
[4478]127DO ilay=1,nlay
128   DO igrid=1,ngrid
129      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
130   ENDDO
131ENDDO
[4449]132
[4478]133
[4449]134! Computes the gradient Richardson's number and stability functions:
135!===================================================================
136
[4481]137! calculation of cn = Sm value at Ri=0
138! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
139cn=(1./sqrt(cepsilon))**(2/3)
140! calculation of Ri0 such that continuity in slope of Sm at Ri=0
141Ri0=2./rpi*(cinf - cn)*ric/cn
142! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
143Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope
144
145
[4478]146DO ilay=2,nlay
[4449]147    DO igrid=1,ngrid
[4478]148        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
149        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
[4631]150        N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
151        shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
152            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 )
153        Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10)
[4481]154       
155        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
[4478]156            Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
157            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut
[4481]158        ELSE ! stable cases
[4644]159            Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric))
[4478]160            Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope
[4481]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
[4449]165        END IF
166       
167        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
168
169    ENDDO
170ENDDO
171
[4631]172
173! Computing the mixing length:
174!==============================================================
175
176switch_num(:,:)=.false.
177
178IF (iflag_atke_lmix .EQ. 1 ) THEN
179
180   DO ilay=2,nlay
181      DO igrid=1,ngrid
182          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
183          IF (N2(igrid,ilay) .GT. 0.) THEN
184             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
185             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
186                l_exchange(igrid,ilay)=max(lstrat,lmin)
187                switch_num(igrid,ilay)=.true.
188             ENDIF
189          ENDIF
190      ENDDO
191   ENDDO
192
[4632]193ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
194! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms, eq 2
195DO ilay=2,nlay
196      DO igrid=1,ngrid
197          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
198          IF (N2(igrid,ilay) .GT. 0.) THEN
[4644]199             lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1.E-10)*(1.+sqrt(Ri(igrid,ilay))/2.))
[4632]200             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
201                l_exchange(igrid,ilay)=max(lstrat,lmin)
202             ENDIF
203          ENDIF
204      ENDDO
205   ENDDO
206
207
208
[4631]209ELSE
[4632]210! default: neglect effect of local stratification and shear
[4631]211
212   DO ilay=2,nlay+1
213      DO igrid=1,ngrid
214          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
215      ENDDO
216
217   ENDDO
218ENDIF
219
220
[4644]221! Computing the TKE k>=2:
222!========================
[4449]223IF (iflag_atke == 0) THEN
224
[4644]225! stationary solution (dtke/dt=0)
226
[4631]227   DO ilay=2,nlay
[4449]228        DO igrid=1,ngrid
229            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
[4631]230            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4449]231        ENDDO
232    ENDDO
233
[4631]234ELSE IF (iflag_atke == 1) THEN
235
236! full implicit scheme resolved with a second order polynomial equation
237
238    DO ilay=2,nlay
239        DO igrid=1,ngrid
[4644]240           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
[4631]241           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
242                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) &
243                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
244           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
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
[4653]301ELSE IF (iflag_atke == 5) THEN
302! semi implicit scheme from Arpege (V. Masson methodology with
303! Taylor expansion of the dissipation term)
304! and implicit resolution when switch num (when we use the mixing length as a function of tke)
305    DO ilay=2,nlay
306        DO igrid=1,ngrid
307           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
308           if (switch_num(igrid,ilay) .and. N2(igrid,ilay)>0.) then
309             invtau=clmix*Sm(igrid,ilay)/sqrt(N2(igrid,ilay))*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))  &
310                   -sqrt(N2(igrid,ilay))/(cepsilon*clmix)
311             !qq=qq/(1.-dtime*invtau)
312             !qq=qq*exp(dtime*invtau)
313             tke(igrid,ilay)=tke(igrid,ilay)/(1.-dtime*invtau)
314             tke(igrid,ilay)=max(0.,tke(igrid,ilay))
315           else
316             qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
317             +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
318             /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
319             qq=max(0.,qq)
320             tke(igrid,ilay)=0.5*(qq**2)
321           endif
322        ENDDO
323    ENDDO
324
325
[4631]326ELSE
[4463]327   call abort_physic("atke_compute_km_kh", &
[4631]328        'numerical treatment of TKE not possible yet', 1)
[4449]329
330END IF
331
[4644]332! We impose a 0 tke at nlay+1
333!==============================
[4449]334
[4644]335DO igrid=1,ngrid
336 tke(igrid,nlay+1)=0.
337END DO
338
339
340! Calculation of surface TKE (k=1)
341!=================================
342! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
343DO igrid=1,ngrid
344 ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
345 tke(igrid,1)=ctkes*(ustar**2)
346END DO
347
348
349! vertical diffusion of TKE
350!==========================
[4653]351IF (atke_ok_vdiff) THEN
[4644]352   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
353ENDIF
354
355
[4449]356! Computing eddy diffusivity coefficients:
357!========================================
[4478]358DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
[4449]359    DO igrid=1,ngrid
[4478]360        ! we add the molecular viscosity to Km,h
361        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
362        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
[4449]363    END DO
364END DO
365
[4478]366! for output:
367!===========
368Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
369Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
[4449]370
371end subroutine atke_compute_km_kh
372
[4644]373!===============================================================================================
374subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
[4449]375
[4644]376! routine that computes the vertical diffusion of TKE by the turbulence
[4653]377! using an implicit resolution (See note by Dufresne and Ghattas (2009))
[4644]378! E Vignon, July 2023
379
380USE atke_turbulence_ini_mod, ONLY : rd, cke, viscom
381
382
383INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
384INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
385
386REAL, INTENT(IN)    :: dtime ! physics time step (s)
387REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
388REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
389REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
390REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
391REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: l_exchange     ! mixing length at interfaces between layers
392REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
393
394REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
395
396
397
398INTEGER                                       :: igrid,ilay
399REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
400REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
401REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
402REAL                                          :: gammak,Kem,KKb,KKt
403
404
405! Few initialisations
406CCK(:,:)=0.
407DDK(:,:)=0.
408dtke(:,:)=0.
409
410
411! Eddy diffusivity for TKE
412
413DO ilay=2,nlay
414    DO igrid=1,ngrid
415       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke
416    ENDDO
417ENDDO
418! at the top of the atmosphere set to 0
419Ke(:,nlay+1)=0.
420! at the surface, set it equal to that at the first model level
421Ke(:,1)=Ke(:,2)
422
423
424! calculate intermediary variables
425
426DO ilay=2,nlay
427    DO igrid=1,ngrid
428    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))   
429    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
430    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
431    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
432    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
433    ak(igrid,ilay)=-gammak*dtime*KKb
434    ck(igrid,ilay)=-gammak*dtime*KKt
435    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
436    ENDDO
437ENDDO
438
439! calculate CCK and DDK coefficients
440! downhill phase
441
442DO igrid=1,ngrid
443  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
444  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
445ENDDO
446
447
448DO ilay=nlay-1,2,-1
449    DO igrid=1,ngrid
450        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
451                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
452        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
453    ENDDO
454ENDDO
455
456! calculate TKE
457! uphill phase
458
459DO ilay=2,nlay+1
460    DO igrid=1,ngrid
461        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
462    ENDDO
463ENDDO
464
465! update TKE
466tke(:,:)=tke(:,:)+dtke(:,:)
467
468
469end subroutine atke_vdiff_tke
470
471
472
[4449]473end module atke_exchange_coeff_mod
Note: See TracBrowser for help on using the repository browser.