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

Last change on this file since 4636 was 4632, checked in by evignon, 12 months ago

petits ajustements suite a la derniere commission
File size: 11.1 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 IF (iflag_atke_lmix .EQ. 2 ) THEN
179! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms, eq 2
180DO ilay=2,nlay
181      DO igrid=1,ngrid
182          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
183          IF (N2(igrid,ilay) .GT. 0.) THEN
184             lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1E-10)*(1.+sqrt(Ri(igrid,ilay))/2.))
185             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
186                l_exchange(igrid,ilay)=max(lstrat,lmin)
187             ENDIF
188          ENDIF
189      ENDDO
190   ENDDO
191
192
193
194ELSE
195! default: neglect effect of local stratification and shear
196
197   DO ilay=2,nlay+1
198      DO igrid=1,ngrid
199          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
200      ENDDO
201
202   ENDDO
203ENDIF
204
205
206! Computing the TKE:
207!===================
208IF (iflag_atke == 0) THEN
209
210   DO ilay=2,nlay
211        DO igrid=1,ngrid
212            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
213            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
214        ENDDO
215    ENDDO
216
217ELSE IF (iflag_atke == 1) THEN
218
219! full implicit scheme resolved with a second order polynomial equation
220
221    DO ilay=2,nlay
222        DO igrid=1,ngrid
223           qq=sqrt(2.*tke(igrid,ilay))
224           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
225                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) &
226                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
227           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
228           tke(igrid,ilay)=0.5*(qq**2)
229        ENDDO
230    ENDDO
231
232ELSE IF (iflag_atke == 2) THEN
233
234! full implicit scheme resolved with a second order polynomial equation
235! + exact resolution for very stable cases (iflag_atke_lmix=1)
236    DO ilay=2,nlay
237        DO igrid=1,ngrid
238           qq=sqrt(2.*tke(igrid,ilay))
239           IF (switch_num(igrid,ilay)) THEN
240           taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
241                    - sqrt(N2(igrid,ilay))/2./cepsilon/clmix
242           taustrat=-1./taustrat
243           qq=qq*exp(-dtime/taustrat)
244           ELSE
245           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
246                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) &
247                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
248           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
249           ENDIF
250           tke(igrid,ilay)=0.5*(qq**2)
251        ENDDO
252    ENDDO
253
254
255
256ELSE IF (iflag_atke == 3) THEN
257
258! semi implicit scheme when l does not depend on tke
259! positive-guaranteed if pr slope in stable condition >1
260    DO ilay=2,nlay
261        DO igrid=1,ngrid
262           qq=sqrt(2.*tke(igrid,ilay))
263           qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
264             /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)))     
265           tke(igrid,ilay)=0.5*(qq**2)
266        ENDDO
267    ENDDO
268
269
270ELSE
271   call abort_physic("atke_compute_km_kh", &
272        'numerical treatment of TKE not possible yet', 1)
273
274END IF
275
276
277! Computing eddy diffusivity coefficients:
278!========================================
279DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
280    DO igrid=1,ngrid
281        ! we add the molecular viscosity to Km,h
282        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
283        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
284    END DO
285END DO
286
287! for output:
288!===========
289Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
290Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
291
292end subroutine atke_compute_km_kh
293
294
295end module atke_exchange_coeff_mod
Note: See TracBrowser for help on using the repository browser.