source: LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 @ 4639

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

commission pour l'atelier TKE avec introduction d'une resolution pronostique
de lequation de la tke
File size: 5.3 KB
Line 
1module call_atke_mod
2
3USE atke_exchange_coeff_mod, ONLY :  atke_compute_km_kh
4
5implicit none
6
7
8contains
9
10subroutine 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
17USE atke_turbulence_ini_mod, ONLY : iflag_num_atke, rg, rd
18
19implicit none
20
21
22! Declarations:
23!=============
24
25REAL, INTENT(IN)    :: dtime ! timestep (s)
26INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
27INTEGER, INTENT(IN) :: nlay ! number of vertical index 
28
29
30REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_uv  ! drag coefficient for wind
31REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_t   ! drag coefficient for temperature
32
33REAL, DIMENSION(ngrid), INTENT(IN)       :: u_surf    ! x wind velocity at the surface
34REAL, DIMENSION(ngrid), INTENT(IN)       :: v_surf    ! y wind velocity at the surface
35REAL, DIMENSION(ngrid), INTENT(IN)       :: temp_surf ! surface temperature
36
37REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
38REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
39REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
40REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
41REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
42
43
44REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
45
46REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
47REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers
48
49
50REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict
51REAL, DIMENSION(ngrid) ::  wind1
52INTEGER i
53
54
55
56call atke_compute_km_kh(ngrid,nlay,dtime,&
57                        wind_u,wind_v,temp,play,pinterf, &
58                        tke,Km_out,Kh_out)
59
60if (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,dtime,&
73                        wind_u_predict,wind_v_predict,temp,play,pinterf, &
74                        tke,Km_out,Kh_out)
75
76end if
77
78end subroutine call_atke
79
80!----------------------------------------------------------------------------------------
81
82subroutine atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,x_in,K_in,x_surf,cdrag,x_predict)
83
84INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
85INTEGER, INTENT(IN) :: nlay ! number of vertical index 
86REAL, INTENT(IN) :: rg,rd,dtime ! gravity, R dry air and timestep
87REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: play      ! pressure middle of layers (Pa)
88REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: temp      ! temperature (K)
89REAL, DIMENSION(ngrid),        INTENT(IN)  :: wind1     ! wind speed first level (m/s)
90REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: pinterf ! pressure at interfaces(Pa)
91REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: x_in    ! variable at the beginning of timestep
92REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: K_in    ! eddy diffusivity coef at the beginning of time step
93REAL, DIMENSION(ngrid),        INTENT(IN)  :: x_surf  ! surface variable
94REAL, DIMENSION(ngrid),        INTENT(IN)  :: cdrag  ! drag coefficient
95REAL, DIMENSION(ngrid,nlay),   INTENT(OUT) :: x_predict  ! variable at the end of time step after explicit prediction
96
97
98integer i,k
99real ml,F1,rho
100real, dimension(ngrid) :: play1,temp1
101real, dimension(ngrid,nlay+1) :: K_big
102
103! computation of K_big
104
105play1(:)=play(:,1)
106temp1(:)=temp(:,1)
107
108! "big K" calculation
109do  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
114enddo
115! speficic treatment for k=nlay
116do 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)
119enddo
120
121
122
123! x_predict calculation for 2<=k<=nlay-1
124do  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
131enddo
132
133! Specific treatment for k=1
134do 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)
138enddo
139
140! Specific treatment for k=nlay
141! flux at the top of the atmosphere=0
142do 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)))
145enddo
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
156end subroutine atke_explicit_prediction
157
158
159
160end module call_atke_mod
Note: See TracBrowser for help on using the repository browser.