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

Last change on this file since 4780 was 4780, checked in by evignon, 5 months ago

petites modifs dans les routines atke pour etre en accord avec les notations du papier
File size: 17.1 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, ri0, ri1
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, cn
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    :: preff      ! reference pressure for potential temperature calculations
72REAL    :: thetam     ! mean potential temperature at interface
73REAL    :: delta      ! discriminant of the second order polynomial
74REAL    :: qq         ! tke=qq**2/2
75REAL    :: shear      ! wind shear
76REAL    :: lstrat     ! mixing length depending on local stratification
77REAL    :: taustrat   ! caracteristic timescale for turbulence in very stable conditions
78REAL    :: netloss    ! net loss term of tke
79REAL    :: netsource  ! net source term of tke
80REAL    :: ustar      ! friction velocity estimation
81REAL    :: invtau     
82REAL    :: rvap
83
84! Initializations:
85!================
86
87DO igrid=1,ngrid
88    dz_interf(igrid,1) = 0.0
89    z_interf(igrid,1) = 0.0
90END DO
91
92! Calculation of potential temperature: (if vapor -> virtual potential temperature)
93!=====================================
94
95preff=100000.
96! results should not depend on the choice of preff
97DO ilay=1,nlay
98     DO igrid = 1, ngrid
99        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
100     END DO
101END DO
102
103! account for water vapor mass for buoyancy calculation
104IF (atke_ok_virtual) THEN
105  DO ilay=1,nlay
106     DO igrid = 1, ngrid
107        rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
108        theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
109     END DO
110  END DO
111ENDIF
112
113
114! Calculation of altitude of layers' middle and bottom interfaces:
115!=================================================================
116
117DO ilay=2,nlay+1
118    DO igrid=1,ngrid
119        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
120        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
121    ENDDO
122ENDDO
123
124DO ilay=1,nlay
125   DO igrid=1,ngrid
126      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
127   ENDDO
128ENDDO
129
130
131! Computes the gradient Richardson's number and stability functions:
132!===================================================================
133
134DO ilay=2,nlay
135    DO igrid=1,ngrid
136        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
137        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
138        N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
139        shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
140            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 )
141        Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10)
142       
143        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
144            Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
145            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut
146        ELSE ! stable cases
147            Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric))
148            ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
149            Prandtl(igrid,ilay) = pr_neut*exp(-pr_slope/pr_neut*Ri(igrid,ilay)+Ri(igrid,ilay)/pr_neut) &
150                                + Ri(igrid,ilay) * pr_slope
151            IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
152               call abort_physic("atke_compute_km_kh", &
153               'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
154            ENDIF
155        END IF
156       
157        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
158
159    ENDDO
160ENDDO
161
162
163! Computing the mixing length:
164!==============================================================
165
166
167IF (iflag_atke_lmix .EQ. 1 ) THEN
168! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions
169   DO ilay=2,nlay
170      DO igrid=1,ngrid
171          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
172          IF (N2(igrid,ilay) .GT. 0.) THEN
173             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
174             lstrat=max(lstrat,lmin)
175             !Inverse interpolation, Van de Wiel et al. 2010
176             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
177          ENDIF
178      ENDDO
179   ENDDO
180
181ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
182! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
183! implies 2 tuning coefficients clmix and clmixshear
184DO ilay=2,nlay
185      DO igrid=1,ngrid
186          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
187          IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
188             lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &
189                    clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay)))
190             lstrat=max(lstrat,lmin)
191             !Inverse interpolation, Van de Wiel et al. 2010   
192             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
193          ENDIF
194      ENDDO
195   ENDDO
196
197ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN
198! add effect of wind shear on lstrat following grisogono 2010, qjrms
199! keeping a single tuning coefficient clmix
200DO ilay=2,nlay
201      DO igrid=1,ngrid
202          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
203          IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
204             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay)))
205             lstrat=max(lstrat,lmin)
206             !Inverse interpolation, Van de Wiel et al. 2010   
207             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
208          ENDIF
209      ENDDO
210   ENDDO
211
212
213
214ELSE
215! default Blackadar formulation: neglect effect of local stratification and shear
216
217   DO ilay=2,nlay+1
218      DO igrid=1,ngrid
219          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
220      ENDDO
221
222   ENDDO
223ENDIF
224
225
226! Computing the TKE k>=2:
227!========================
228IF (iflag_atke == 0) THEN
229
230! stationary solution (dtke/dt=0)
231
232   DO ilay=2,nlay
233        DO igrid=1,ngrid
234            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
235            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
236        ENDDO
237    ENDDO
238
239ELSE IF (iflag_atke == 1) THEN
240
241! full implicit scheme resolved with a second order polynomial equation
242! default solution which shows fair convergence properties
243    DO ilay=2,nlay
244        DO igrid=1,ngrid
245           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
246           delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. &
247                 +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + &
248                 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) &
249                 *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
250           qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
251           qq=max(0.,qq)
252           tke(igrid,ilay)=0.5*(qq**2)
253        ENDDO
254    ENDDO
255
256
257ELSE IF (iflag_atke == 2) THEN
258
259! semi implicit scheme when l does not depend on tke
260! positive-guaranteed if pr slope in stable condition >1
261
262   DO ilay=2,nlay
263        DO igrid=1,ngrid
264           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
265           qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      &
266               *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
267               /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
268           tke(igrid,ilay)=0.5*(qq**2)
269        ENDDO
270    ENDDO
271
272
273ELSE IF (iflag_atke == 3) THEN
274! numerical resolution adapted from that in MAR (Deleersnijder 1992)
275! positively defined by construction
276
277    DO ilay=2,nlay
278        DO igrid=1,ngrid
279           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
280           IF (Ri(igrid,ilay) .LT. 0.) THEN
281              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
282              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
283           ELSE
284              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
285                      l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
286              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
287           ENDIF
288           qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
289           tke(igrid,ilay)=0.5*(qq**2)
290        ENDDO
291    ENDDO
292
293ELSE IF (iflag_atke == 4) THEN
294! semi implicit scheme from Arpege (V. Masson methodology with
295! Taylor expansion of the dissipation term)
296    DO ilay=2,nlay
297        DO igrid=1,ngrid
298           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
299           qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
300             +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
301             /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
302           qq=max(0.,qq)
303           tke(igrid,ilay)=0.5*(qq**2)
304        ENDDO
305    ENDDO
306
307
308ELSE
309   call abort_physic("atke_compute_km_kh", &
310        'numerical treatment of TKE not possible yet', 1)
311
312END IF
313
314! We impose a 0 tke at nlay+1
315!==============================
316
317DO igrid=1,ngrid
318 tke(igrid,nlay+1)=0.
319END DO
320
321
322! Calculation of surface TKE (k=1)
323!=================================
324! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
325DO igrid=1,ngrid
326 ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
327 tke(igrid,1)=ctkes*(ustar**2)
328END DO
329
330
331! vertical diffusion of TKE
332!==========================
333IF (atke_ok_vdiff) THEN
334   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
335ENDIF
336
337
338! Computing eddy diffusivity coefficients:
339!========================================
340DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
341    DO igrid=1,ngrid
342        ! we add the molecular viscosity to Km,h
343        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
344        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
345    END DO
346END DO
347
348! for output:
349!===========
350Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
351Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
352
353end subroutine atke_compute_km_kh
354
355!===============================================================================================
356subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
357
358! routine that computes the vertical diffusion of TKE by the turbulence
359! using an implicit resolution (See note by Dufresne and Ghattas (2009))
360! E Vignon, July 2023
361
362USE lmdz_atke_turbulence_ini, ONLY : rd, cke, viscom
363
364
365INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
366INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
367
368REAL, INTENT(IN)    :: dtime ! physics time step (s)
369REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
370REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
371REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
372REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
373REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: l_exchange     ! mixing length at interfaces between layers
374REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
375
376REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
377
378
379
380INTEGER                                       :: igrid,ilay
381REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
382REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
383REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
384REAL                                          :: gammak,Kem,KKb,KKt
385
386
387! Few initialisations
388CCK(:,:)=0.
389DDK(:,:)=0.
390dtke(:,:)=0.
391
392
393! Eddy diffusivity for TKE
394
395DO ilay=2,nlay
396    DO igrid=1,ngrid
397       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke
398    ENDDO
399ENDDO
400! at the top of the atmosphere set to 0
401Ke(:,nlay+1)=0.
402! at the surface, set it equal to that at the first model level
403Ke(:,1)=Ke(:,2)
404
405
406! calculate intermediary variables
407
408DO ilay=2,nlay
409    DO igrid=1,ngrid
410    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))   
411    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
412    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
413    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
414    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
415    ak(igrid,ilay)=-gammak*dtime*KKb
416    ck(igrid,ilay)=-gammak*dtime*KKt
417    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
418    ENDDO
419ENDDO
420
421! calculate CCK and DDK coefficients
422! downhill phase
423
424DO igrid=1,ngrid
425  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
426  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
427ENDDO
428
429
430DO ilay=nlay-1,2,-1
431    DO igrid=1,ngrid
432        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
433                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
434        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
435    ENDDO
436ENDDO
437
438! calculate TKE
439! uphill phase
440
441DO ilay=2,nlay+1
442    DO igrid=1,ngrid
443        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
444    ENDDO
445ENDDO
446
447! update TKE
448tke(:,:)=tke(:,:)+dtke(:,:)
449
450
451end subroutine atke_vdiff_tke
452
453
454
455end module lmdz_atke_exchange_coeff
Note: See TracBrowser for help on using the repository browser.