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

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

commission pour l'atelier TKE avec introduction d'une resolution pronostique
de lequation de la tke
File size: 10.5 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,play,pinterf, &
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
29USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd
30USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin
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)       :: play   ! pressure (Pa)
46REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
47
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)
69LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases
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
80
81! Initializations:
82!================
83
84DO igrid=1,ngrid
85    dz_interf(igrid,1) = 0.0
86    z_interf(igrid,1) = 0.0
87END DO
88
89! Calculation of potential temperature: (if vapor -> todo virtual potential temperature)
90!=====================================
91
92preff=100000.
93! The result should not depend on the choice of preff
94DO ilay=1,nlay
95     DO igrid = 1, ngrid
96        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
97     END DO
98END DO
99
100
101
102! Calculation of altitude of layers' middle and bottom interfaces:
103!=================================================================
104
105DO ilay=2,nlay+1
106    DO igrid=1,ngrid
107        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
108        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
109    ENDDO
110ENDDO
111
112DO ilay=1,nlay
113   DO igrid=1,ngrid
114      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
115   ENDDO
116ENDDO
117
118
119! Computes the gradient Richardson's number and stability functions:
120!===================================================================
121
122! calculation of cn = Sm value at Ri=0
123! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
124cn=(1./sqrt(cepsilon))**(2/3)
125! calculation of Ri0 such that continuity in slope of Sm at Ri=0
126Ri0=2./rpi*(cinf - cn)*ric/cn
127! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
128Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope
129
130
131DO ilay=2,nlay
132    DO igrid=1,ngrid
133        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
134        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
135        N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
136        shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
137            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 )
138        Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10)
139       
140        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
141            Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
142            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut
143        ELSE ! stable cases
144            Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric))
145            Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope
146            IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
147               call abort_physic("atke_compute_km_kh", &
148               'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
149            ENDIF
150        END IF
151       
152        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
153
154    ENDDO
155ENDDO
156
157
158! Computing the mixing length:
159!==============================================================
160
161switch_num(:,:)=.false.
162
163IF (iflag_atke_lmix .EQ. 1 ) THEN
164
165   DO ilay=2,nlay
166      DO igrid=1,ngrid
167          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
168          IF (N2(igrid,ilay) .GT. 0.) THEN
169             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
170             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
171                l_exchange(igrid,ilay)=max(lstrat,lmin)
172                switch_num(igrid,ilay)=.true.
173             ENDIF
174          ENDIF
175      ENDDO
176   ENDDO
177
178ELSE
179! default: neglect effect of local stratification
180
181   DO ilay=2,nlay+1
182      DO igrid=1,ngrid
183          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
184      ENDDO
185
186   ENDDO
187ENDIF
188
189
190! Computing the TKE:
191!===================
192IF (iflag_atke == 0) THEN
193
194   DO ilay=2,nlay
195        DO igrid=1,ngrid
196            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
197            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
198        ENDDO
199    ENDDO
200
201ELSE IF (iflag_atke == 1) THEN
202
203! full implicit scheme resolved with a second order polynomial equation
204
205    DO ilay=2,nlay
206        DO igrid=1,ngrid
207           qq=sqrt(2.*tke(igrid,ilay))
208           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
209                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) &
210                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
211           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
212           tke(igrid,ilay)=0.5*(qq**2)
213        ENDDO
214    ENDDO
215
216ELSE IF (iflag_atke == 2) THEN
217
218! full implicit scheme resolved with a second order polynomial equation
219! + exact resolution for very stable cases
220    DO ilay=2,nlay
221        DO igrid=1,ngrid
222           qq=sqrt(2.*tke(igrid,ilay))
223           IF (switch_num(igrid,ilay)) THEN
224           taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
225                    - sqrt(N2(igrid,ilay))/2./cepsilon/clmix
226           taustrat=-1./taustrat
227           qq=qq*exp(-dtime/taustrat)
228           ELSE
229           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
230                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) &
231                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
232           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
233           ENDIF
234           tke(igrid,ilay)=0.5*(qq**2)
235        ENDDO
236    ENDDO
237
238
239
240ELSE IF (iflag_atke == 3) THEN
241
242! semi implicit scheme when l does not depend on tke
243! positive-guaranteed if pr slope in stable condition >1
244    DO ilay=2,nlay
245        DO igrid=1,ngrid
246           qq=sqrt(2.*tke(igrid,ilay))
247           qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
248             /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)))     
249           tke(igrid,ilay)=0.5*(qq**2)
250        ENDDO
251    ENDDO
252
253
254ELSE
255   call abort_physic("atke_compute_km_kh", &
256        'numerical treatment of TKE not possible yet', 1)
257
258END IF
259
260
261! Computing eddy diffusivity coefficients:
262!========================================
263DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
264    DO igrid=1,ngrid
265        ! we add the molecular viscosity to Km,h
266        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
267        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
268    END DO
269END DO
270
271! for output:
272!===========
273Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
274Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
275
276end subroutine atke_compute_km_kh
277
278
279end module atke_exchange_coeff_mod
Note: See TracBrowser for help on using the repository browser.