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

Last change on this file since 4799 was 4780, checked in by evignon, 13 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.