source: trunk/LMDZ.MARS/libf/phymars/lmdz_atke_exchange_coeff.F90 @ 3567

Last change on this file since 3567 was 3168, checked in by llange, 12 months ago

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