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

Last change on this file since 4678 was 4675, checked in by evignon, 9 months ago

correction option numerique 1 pour ATKE
File size: 16.7 KB
Line 
1module atke_exchange_coeff_mod
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) dry atmosphere
23! (2) horizontal homogeneity (Dx=Dy=0.)
24!=======================================================================
25
26
27
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
30USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin
31
32implicit none
33
34
35! Declarations:
36!=============
37
38INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
39INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
40
41REAL, INTENT(IN)    :: dtime ! physics time step (s)
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)
45REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap   ! specific humidity (kg/kg)
46REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
47REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
48REAL, DIMENSION(ngrid), INTENT(IN)            :: cdrag_uv   ! surface drag coefficient for momentum
49
50REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
51
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
54
55! Local variables
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)
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)
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)
70
71INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
72REAL    :: cn,Ri0,Ri1    ! parameter for Sm stability function and Prandlt
73REAL    :: preff      ! reference pressure for potential temperature calculations
74REAL    :: thetam     ! mean potential temperature at interface
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
80REAL    :: netloss    ! net loss term of tke
81REAL    :: netsource  ! net source term of tke
82REAL    :: ustar      ! friction velocity estimation
83REAL    :: invtau     
84REAL    :: rvap
85
86! Initializations:
87!================
88
89DO igrid=1,ngrid
90    dz_interf(igrid,1) = 0.0
91    z_interf(igrid,1) = 0.0
92END DO
93
94! Calculation of potential temperature: (if vapor -> virtual potential temperature)
95!=====================================
96
97preff=100000.
98! results should not depend on the choice of preff
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
104
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
114
115
116! Calculation of altitude of layers' middle and bottom interfaces:
117!=================================================================
118
119DO ilay=2,nlay+1
120    DO igrid=1,ngrid
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)
123    ENDDO
124ENDDO
125
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
131
132
133! Computes the gradient Richardson's number and stability functions:
134!===================================================================
135
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
145DO ilay=2,nlay
146    DO igrid=1,ngrid
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))
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)
153       
154        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
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
157        ELSE ! stable cases
158            Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric))
159            Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope
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
164        END IF
165       
166        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
167
168    ENDDO
169ENDDO
170
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))
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)
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             !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
204          ENDIF
205      ENDDO
206   ENDDO
207
208
209
210ELSE
211! default: neglect effect of local stratification and shear
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
222! Computing the TKE k>=2:
223!========================
224IF (iflag_atke == 0) THEN
225
226! stationary solution (dtke/dt=0)
227
228   DO ilay=2,nlay
229        DO igrid=1,ngrid
230            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
231            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
232        ENDDO
233    ENDDO
234
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
241           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
242           delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. &
243                 +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + &
244                 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) &
245                 *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
246           qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
247           qq=max(0.,qq)
248           tke(igrid,ilay)=0.5*(qq**2)
249        ENDDO
250    ENDDO
251
252
253ELSE IF (iflag_atke == 2) THEN
254
255! semi implicit scheme when l does not depend on tke
256! positive-guaranteed if pr slope in stable condition >1
257
258   DO ilay=2,nlay
259        DO igrid=1,ngrid
260           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
261           qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      &
262               *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
263               /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
264           tke(igrid,ilay)=0.5*(qq**2)
265        ENDDO
266    ENDDO
267
268
269ELSE IF (iflag_atke == 3) THEN
270! numerical resolution adapted from that in MAR (Deleersnijder 1992)
271! positively defined by construction
272
273    DO ilay=2,nlay
274        DO igrid=1,ngrid
275           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
276           IF (Ri(igrid,ilay) .LT. 0.) THEN
277              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
278              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
279           ELSE
280              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
281                      l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
282              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
283           ENDIF
284           qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
285           tke(igrid,ilay)=0.5*(qq**2)
286        ENDDO
287    ENDDO
288
289ELSE IF (iflag_atke == 4) THEN
290! semi implicit scheme from Arpege (V. Masson methodology with
291! Taylor expansion of the dissipation term)
292    DO ilay=2,nlay
293        DO igrid=1,ngrid
294           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
295           qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
296             +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
297             /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
298           qq=max(0.,qq)
299           tke(igrid,ilay)=0.5*(qq**2)
300        ENDDO
301    ENDDO
302
303
304ELSE
305   call abort_physic("atke_compute_km_kh", &
306        'numerical treatment of TKE not possible yet', 1)
307
308END IF
309
310! We impose a 0 tke at nlay+1
311!==============================
312
313DO igrid=1,ngrid
314 tke(igrid,nlay+1)=0.
315END DO
316
317
318! Calculation of surface TKE (k=1)
319!=================================
320! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
321DO igrid=1,ngrid
322 ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
323 tke(igrid,1)=ctkes*(ustar**2)
324END DO
325
326
327! vertical diffusion of TKE
328!==========================
329IF (atke_ok_vdiff) THEN
330   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
331ENDIF
332
333
334! Computing eddy diffusivity coefficients:
335!========================================
336DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
337    DO igrid=1,ngrid
338        ! we add the molecular viscosity to Km,h
339        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
340        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
341    END DO
342END DO
343
344! for output:
345!===========
346Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
347Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
348
349end subroutine atke_compute_km_kh
350
351!===============================================================================================
352subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
353
354! routine that computes the vertical diffusion of TKE by the turbulence
355! using an implicit resolution (See note by Dufresne and Ghattas (2009))
356! E Vignon, July 2023
357
358USE atke_turbulence_ini_mod, ONLY : rd, cke, viscom
359
360
361INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
362INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
363
364REAL, INTENT(IN)    :: dtime ! physics time step (s)
365REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
366REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
367REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
368REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
369REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: l_exchange     ! mixing length at interfaces between layers
370REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
371
372REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
373
374
375
376INTEGER                                       :: igrid,ilay
377REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
378REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
379REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
380REAL                                          :: gammak,Kem,KKb,KKt
381
382
383! Few initialisations
384CCK(:,:)=0.
385DDK(:,:)=0.
386dtke(:,:)=0.
387
388
389! Eddy diffusivity for TKE
390
391DO ilay=2,nlay
392    DO igrid=1,ngrid
393       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke
394    ENDDO
395ENDDO
396! at the top of the atmosphere set to 0
397Ke(:,nlay+1)=0.
398! at the surface, set it equal to that at the first model level
399Ke(:,1)=Ke(:,2)
400
401
402! calculate intermediary variables
403
404DO ilay=2,nlay
405    DO igrid=1,ngrid
406    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))   
407    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
408    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
409    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
410    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
411    ak(igrid,ilay)=-gammak*dtime*KKb
412    ck(igrid,ilay)=-gammak*dtime*KKt
413    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
414    ENDDO
415ENDDO
416
417! calculate CCK and DDK coefficients
418! downhill phase
419
420DO igrid=1,ngrid
421  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
422  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
423ENDDO
424
425
426DO ilay=nlay-1,2,-1
427    DO igrid=1,ngrid
428        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
429                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
430        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
431    ENDDO
432ENDDO
433
434! calculate TKE
435! uphill phase
436
437DO ilay=2,nlay+1
438    DO igrid=1,ngrid
439        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
440    ENDDO
441ENDDO
442
443! update TKE
444tke(:,:)=tke(:,:)+dtke(:,:)
445
446
447end subroutine atke_vdiff_tke
448
449
450
451end module atke_exchange_coeff_mod
Note: See TracBrowser for help on using the repository browser.