1 |
|
|
module call_atke_mod |
2 |
|
|
|
3 |
|
|
USE atke_exchange_coeff_mod, ONLY : atke_compute_km_kh |
4 |
|
|
|
5 |
|
|
implicit none |
6 |
|
|
|
7 |
|
|
|
8 |
|
|
contains |
9 |
|
|
|
10 |
|
|
subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, & |
11 |
|
|
wind_u,wind_v,temp,play,pinterf, & |
12 |
|
|
tke,Km_out,Kh_out) |
13 |
|
|
|
14 |
|
|
|
15 |
|
|
|
16 |
|
|
|
17 |
|
|
USE atke_turbulence_ini_mod, ONLY : iflag_num_atke, rg, rd |
18 |
|
|
|
19 |
|
|
implicit none |
20 |
|
|
|
21 |
|
|
|
22 |
|
|
! Declarations: |
23 |
|
|
!============= |
24 |
|
|
|
25 |
|
|
REAL, INTENT(IN) :: dtime ! timestep (s) |
26 |
|
|
INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
27 |
|
|
INTEGER, INTENT(IN) :: nlay ! number of vertical index |
28 |
|
|
|
29 |
|
|
|
30 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_uv ! drag coefficient for wind |
31 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_t ! drag coefficient for temperature |
32 |
|
|
|
33 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: u_surf ! x wind velocity at the surface |
34 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: v_surf ! y wind velocity at the surface |
35 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: temp_surf ! surface temperature |
36 |
|
|
|
37 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s) |
38 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s) |
39 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
40 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
41 |
|
|
REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) |
42 |
|
|
|
43 |
|
|
|
44 |
|
|
REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
45 |
|
|
|
46 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers |
47 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers |
48 |
|
|
|
49 |
|
|
|
50 |
|
|
REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict |
51 |
|
|
REAL, DIMENSION(ngrid) :: wind1 |
52 |
|
|
INTEGER i |
53 |
|
|
|
54 |
|
|
|
55 |
|
|
|
56 |
|
|
call atke_compute_km_kh(ngrid,nlay, & |
57 |
|
|
wind_u,wind_v,temp,play,pinterf, & |
58 |
|
|
tke,Km_out,Kh_out) |
59 |
|
|
|
60 |
|
|
if (iflag_num_atke .EQ. 1) then |
61 |
|
|
|
62 |
|
|
!! pay attention that the treatment of the TKE |
63 |
|
|
!! has to be adapted when solving the TKE with a prognostic equation |
64 |
|
|
|
65 |
|
|
do i=1,ngrid |
66 |
|
|
wind1(i)=sqrt(wind_u(i,1)**2+wind_v(i,1)**2) |
67 |
|
|
enddo |
68 |
|
|
call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_u,Km_out,u_surf,cdrag_uv,wind_u_predict) |
69 |
|
|
call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_v,Km_out,v_surf,cdrag_uv,wind_v_predict) |
70 |
|
|
|
71 |
|
|
|
72 |
|
|
call atke_compute_km_kh(ngrid,nlay, & |
73 |
|
|
wind_u_predict,wind_v_predict,temp,play,pinterf, & |
74 |
|
|
tke,Km_out,Kh_out) |
75 |
|
|
|
76 |
|
|
end if |
77 |
|
|
|
78 |
|
|
end subroutine call_atke |
79 |
|
|
|
80 |
|
|
!---------------------------------------------------------------------------------------- |
81 |
|
|
|
82 |
|
|
subroutine atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,x_in,K_in,x_surf,cdrag,x_predict) |
83 |
|
|
|
84 |
|
|
INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
85 |
|
|
INTEGER, INTENT(IN) :: nlay ! number of vertical index |
86 |
|
|
REAL, INTENT(IN) :: rg,rd,dtime ! gravity, R dry air and timestep |
87 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure middle of layers (Pa) |
88 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
89 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: wind1 ! wind speed first level (m/s) |
90 |
|
|
REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) |
91 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: x_in ! variable at the beginning of timestep |
92 |
|
|
REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: K_in ! eddy diffusivity coef at the beginning of time step |
93 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: x_surf ! surface variable |
94 |
|
|
REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag ! drag coefficient |
95 |
|
|
REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: x_predict ! variable at the end of time step after explicit prediction |
96 |
|
|
|
97 |
|
|
|
98 |
|
|
integer i,k |
99 |
|
|
real ml,F1,rho |
100 |
|
|
real, dimension(ngrid) :: play1,temp1 |
101 |
|
|
real, dimension(ngrid,nlay+1) :: K_big |
102 |
|
|
|
103 |
|
|
! computation of K_big |
104 |
|
|
|
105 |
|
|
play1(:)=play(:,1) |
106 |
|
|
temp1(:)=temp(:,1) |
107 |
|
|
|
108 |
|
|
! "big K" calculation |
109 |
|
|
do k=2,nlay-1 |
110 |
|
|
do i=1,ngrid |
111 |
|
|
rho=pinterf(i,k)/rd/(0.5*(temp(i,k-1)+temp(i,k))) |
112 |
|
|
K_big(i,k)=rg*K_in(i,k)/(play(i,k)-play(i,k+1))*(rho**2) |
113 |
|
|
enddo |
114 |
|
|
enddo |
115 |
|
|
! speficic treatment for k=nlay |
116 |
|
|
do i=1,ngrid |
117 |
|
|
rho=pinterf(i,nlay)/rd/temp(i,nlay) |
118 |
|
|
K_big(i,nlay)=rg*K_in(i,nlay)/(2*(play(i,nlay)-pinterf(i,nlay+1)))*(rho**2) |
119 |
|
|
enddo |
120 |
|
|
|
121 |
|
|
|
122 |
|
|
|
123 |
|
|
! x_predict calculation for 2<=k<=nlay-1 |
124 |
|
|
do k=2,nlay-1 |
125 |
|
|
do i=1,ngrid |
126 |
|
|
ml=(pinterf(i,k)-pinterf(i,k+1))/rg |
127 |
|
|
x_predict(i,k)=x_in(i,k)-dtime/ml*(-K_big(i,k+1)*x_in(i,k+1) & |
128 |
|
|
+ (K_big(i,k)+K_big(i,k+1))*x_in(i,k) & |
129 |
|
|
- K_big(i,k)*x_in(i,k-1)) |
130 |
|
|
enddo |
131 |
|
|
enddo |
132 |
|
|
|
133 |
|
|
! Specific treatment for k=1 |
134 |
|
|
do i=1,ngrid |
135 |
|
|
ml=(pinterf(i,1)-pinterf(i,2))/rg |
136 |
|
|
F1=-play1(i)/rd/temp1(i)*wind1(i)*cdrag(i)*(x_in(i,1)-x_surf(i)) ! attention convention sens du flux |
137 |
|
|
x_predict(i,1)=x_in(i,1)-dtime/ml*(-K_big(i,2)*(x_in(i,2) - x_in(i,1))-F1) |
138 |
|
|
enddo |
139 |
|
|
|
140 |
|
|
! Specific treatment for k=nlay |
141 |
|
|
! flux at the top of the atmosphere=0 |
142 |
|
|
do i=1,ngrid |
143 |
|
|
ml=0.5*(pinterf(i,nlay)-pinterf(i,nlay+1))/rg |
144 |
|
|
x_predict(i,nlay)=x_in(i,nlay)+dtime/ml*(K_big(i,nlay)*(x_in(i,nlay)-x_in(i,nlay-1))) |
145 |
|
|
enddo |
146 |
|
|
|
147 |
|
|
!K_big(:,1)=0. |
148 |
|
|
!do k=1,nlay |
149 |
|
|
! do i=1,ngrid |
150 |
|
|
! print*, 'youhou', k, x_in(i,k), x_predict(i,k), K_big(i,k) |
151 |
|
|
! end do |
152 |
|
|
!enddo |
153 |
|
|
|
154 |
|
|
|
155 |
|
|
|
156 |
|
|
end subroutine atke_explicit_prediction |
157 |
|
|
|
158 |
|
|
|
159 |
|
|
|
160 |
|
|
end module call_atke_mod |