1 |
|
|
module atke_exchange_coeff_mod |
2 |
|
|
|
3 |
|
|
implicit none |
4 |
|
|
|
5 |
|
|
contains |
6 |
|
|
|
7 |
|
|
subroutine atke_compute_km_kh(ngrid,nlay, & |
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 |
|
|
|
28 |
|
|
USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd |
29 |
|
|
USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd |
30 |
|
|
USE atke_turbulence_ini_mod, ONLY : viscom, viscoh |
31 |
|
|
|
32 |
|
|
implicit none |
33 |
|
|
|
34 |
|
|
|
35 |
|
|
! Declarations: |
36 |
|
|
!============= |
37 |
|
|
|
38 |
|
|
INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
39 |
|
|
INTEGER, INTENT(IN) :: nlay ! number of vertical index |
40 |
|
|
|
41 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s) |
42 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s) |
43 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
44 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
45 |
|
|
REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) |
46 |
|
|
|
47 |
|
|
|
48 |
|
|
REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
49 |
|
|
|
50 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers |
51 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers |
52 |
|
|
|
53 |
|
|
! Local variables |
54 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: Km ! Exchange coefficient for momentum at interface between layers |
55 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: Kh ! Exchange coefficient for heat flux at interface between layers |
56 |
|
|
REAL, DIMENSION(ngrid,nlay) :: theta ! Potential temperature |
57 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange (at interface) |
58 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface |
59 |
|
|
REAL, DIMENSION(ngrid,nlay) :: z_lay ! Altitude of layers |
60 |
|
|
REAL, DIMENSION(ngrid,nlay) :: dz_interf ! distance between two consecutive interfaces |
61 |
|
|
REAL, DIMENSION(ngrid,nlay) :: dz_lay ! distance between two layer middles (NB: first and last are half layers) |
62 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number (at interface) |
63 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number (at interface) |
64 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface) |
65 |
|
|
REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface) |
66 |
|
|
|
67 |
|
|
INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) |
68 |
|
|
REAL :: cn,Ri0,Ri1 ! parameter for Sm stability function and Prandlt |
69 |
|
|
REAL :: preff ! reference pressure for potential temperature calculations |
70 |
|
|
REAL :: thetam ! mean potential temperature at interface |
71 |
|
|
|
72 |
|
|
|
73 |
|
|
! Initializations: |
74 |
|
|
!================ |
75 |
|
|
|
76 |
|
|
DO igrid=1,ngrid |
77 |
|
|
dz_interf(igrid,1) = 0.0 |
78 |
|
|
z_interf(igrid,1) = 0.0 |
79 |
|
|
END DO |
80 |
|
|
|
81 |
|
|
! Calculation of potential temperature: (if vapor -> todo virtual potential temperature) |
82 |
|
|
!===================================== |
83 |
|
|
|
84 |
|
|
preff=100000. |
85 |
|
|
! The result should not depend on the choice of preff |
86 |
|
|
DO ilay=1,nlay |
87 |
|
|
DO igrid = 1, ngrid |
88 |
|
|
theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) |
89 |
|
|
END DO |
90 |
|
|
END DO |
91 |
|
|
|
92 |
|
|
|
93 |
|
|
|
94 |
|
|
! Calculation of altitude of layers' middle and bottom interfaces: |
95 |
|
|
!================================================================= |
96 |
|
|
|
97 |
|
|
DO ilay=2,nlay+1 |
98 |
|
|
DO igrid=1,ngrid |
99 |
|
|
dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay)) |
100 |
|
|
z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1) |
101 |
|
|
ENDDO |
102 |
|
|
ENDDO |
103 |
|
|
|
104 |
|
|
DO ilay=1,nlay |
105 |
|
|
DO igrid=1,ngrid |
106 |
|
|
z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) |
107 |
|
|
ENDDO |
108 |
|
|
ENDDO |
109 |
|
|
|
110 |
|
|
|
111 |
|
|
! Computing the mixing length: |
112 |
|
|
! so far, we have neglected the effect of local stratification |
113 |
|
|
!============================================================== |
114 |
|
|
|
115 |
|
|
|
116 |
|
|
DO ilay=2,nlay+1 |
117 |
|
|
DO igrid=1,ngrid |
118 |
|
|
l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
119 |
|
|
ENDDO |
120 |
|
|
ENDDO |
121 |
|
|
|
122 |
|
|
! Computes the gradient Richardson's number and stability functions: |
123 |
|
|
!=================================================================== |
124 |
|
|
|
125 |
|
|
! calculation of cn = Sm value at Ri=0 |
126 |
|
|
! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0 |
127 |
|
|
cn=(1./sqrt(cepsilon))**(2/3) |
128 |
|
|
! calculation of Ri0 such that continuity in slope of Sm at Ri=0 |
129 |
|
|
Ri0=2./rpi*(cinf - cn)*ric/cn |
130 |
|
|
! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 |
131 |
|
|
Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope |
132 |
|
|
|
133 |
|
|
|
134 |
|
|
DO 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 |
|
|
Ri(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) / & |
139 |
|
|
MAX(((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,1E-10) |
141 |
|
|
|
142 |
|
|
IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases |
143 |
|
|
Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn |
144 |
|
|
Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut |
145 |
|
|
ELSE ! stable cases |
146 |
|
|
Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric)) |
147 |
|
|
Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope |
148 |
|
|
IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN |
149 |
|
|
call abort_physic("atke_compute_km_kh", & |
150 |
|
|
'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1) |
151 |
|
|
ENDIF |
152 |
|
|
END IF |
153 |
|
|
|
154 |
|
|
Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay) |
155 |
|
|
|
156 |
|
|
ENDDO |
157 |
|
|
ENDDO |
158 |
|
|
|
159 |
|
|
! Computing the TKE: |
160 |
|
|
!=================== |
161 |
|
|
IF (iflag_atke == 0) THEN |
162 |
|
|
|
163 |
|
|
! stationary solution neglecting the vertical transport of TKE by turbulence |
164 |
|
|
DO ilay=2,nlay |
165 |
|
|
DO igrid=1,ngrid |
166 |
|
|
tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & |
167 |
|
|
(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & |
168 |
|
|
((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) * & |
169 |
|
|
(1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) |
170 |
|
|
ENDDO |
171 |
|
|
ENDDO |
172 |
|
|
|
173 |
|
|
ELSE ! TO DO |
174 |
|
|
|
175 |
|
|
call abort_physic("atke_compute_km_kh", & |
176 |
|
|
'traitement non-stationnaire de la tke pas encore prevu', 1) |
177 |
|
|
|
178 |
|
|
END IF |
179 |
|
|
|
180 |
|
|
|
181 |
|
|
! Computing eddy diffusivity coefficients: |
182 |
|
|
!======================================== |
183 |
|
|
DO ilay=2,nlay ! TODO: also calculate for nlay+1 ? |
184 |
|
|
DO igrid=1,ngrid |
185 |
|
|
! we add the molecular viscosity to Km,h |
186 |
|
|
Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 |
187 |
|
|
Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 |
188 |
|
|
END DO |
189 |
|
|
END DO |
190 |
|
|
|
191 |
|
|
! for output: |
192 |
|
|
!=========== |
193 |
|
|
Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay) |
194 |
|
|
Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay) |
195 |
|
|
|
196 |
|
|
end subroutine atke_compute_km_kh |
197 |
|
|
|
198 |
|
|
|
199 |
|
|
end module atke_exchange_coeff_mod |