source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_call_atke.F90 @ 5088

Last change on this file since 5088 was 5082, checked in by abarral, 4 months ago

(lint) Fix obsolete boolean operators
File size: 6.2 KB
Line 
1module lmdz_call_atke
2
3USE lmdz_atke_exchange_coeff, ONLY :  atke_compute_km_kh
4
5implicit none
6
7
8contains
9
10subroutine call_atke(dtime,ngrid,nlay,nsrf,ni,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
11                        wind_u,wind_v,temp,qvap,play,pinterf, &
12                        tke,eps,Km_out,Kh_out)
13
14
15
16
17USE lmdz_atke_turbulence_ini, ONLY : iflag_num_atke, rg, rd
18USE phys_local_var_mod, ONLY: tke_shear, tke_buoy, tke_trans
19
20implicit none
21
22
23! Declarations:
24!=============
25
26REAL, INTENT(IN)    :: dtime ! timestep (s)
27INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
28INTEGER, INTENT(IN) :: nlay ! number of vertical index 
29INTEGER, INTENT(IN) :: nsrf ! surface tile index
30INTEGER, DIMENSION(ngrid), INTENT(IN) :: ni ! array of indices to move from knon to klon arrays
31
32
33REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_uv  ! drag coefficient for wind
34REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_t   ! drag coefficient for temperature
35
36REAL, DIMENSION(ngrid), INTENT(IN)       :: u_surf    ! x wind velocity at the surface
37REAL, DIMENSION(ngrid), INTENT(IN)       :: v_surf    ! y wind velocity at the surface
38REAL, DIMENSION(ngrid), INTENT(IN)       :: temp_surf ! surface temperature
39
40REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
41REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
42REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
43REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap   ! specific humidity (kg/kg)
44REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
45REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
46
47
48REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke      ! turbulent kinetic energy at interface between layers
49
50REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: eps      ! output: tke dissipation rate at interface between layers
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
55REAL, DIMENSION(ngrid,nlay+1) :: tke_shear_term,tke_buoy_term,tke_trans_term
56REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict
57REAL, DIMENSION(ngrid) ::  wind1
58INTEGER i,j,k
59
60
61call atke_compute_km_kh(ngrid,nlay,dtime,&
62                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
63                        tke,eps,tke_shear_term,tke_buoy_term,tke_trans_term,Km_out,Kh_out)
64
65
66                               
67if (iflag_num_atke == 1) then
68   !! In this case, we make an explicit prediction of the wind shear to calculate the tke in a
69   !! forward backward way
70   !! pay attention that the treatment of the TKE
71   !! has to be adapted when solving the TKE with a prognostic equation
72
73   do i=1,ngrid
74      wind1(i)=sqrt(wind_u(i,1)**2+wind_v(i,1)**2)
75   enddo
76   call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_u,Km_out,u_surf,cdrag_uv,wind_u_predict) 
77   call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_v,Km_out,v_surf,cdrag_uv,wind_v_predict)
78
79
80   call atke_compute_km_kh(ngrid,nlay,dtime,&
81                        wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
82                        tke,eps,tke_shear_term,tke_buoy_term,tke_trans_term,Km_out,Kh_out)
83
84end if
85
86
87! Diagnostics of tke loss/source terms
88
89 DO k=1,nlay+1
90    DO i=1,ngrid
91       j=ni(i)
92       tke_shear(j,k,nsrf)=tke_shear_term(i,k)
93       tke_buoy(j,k,nsrf)=tke_buoy_term(i,k)
94       tke_trans(j,k,nsrf)=tke_trans_term(i,k)
95    ENDDO
96 ENDDO
97
98
99
100end subroutine call_atke
101
102!----------------------------------------------------------------------------------------
103
104subroutine atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,x_in,K_in,x_surf,cdrag,x_predict)
105
106INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
107INTEGER, INTENT(IN) :: nlay ! number of vertical index 
108REAL, INTENT(IN) :: rg,rd,dtime ! gravity, R dry air and timestep
109REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: play      ! pressure middle of layers (Pa)
110REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: temp      ! temperature (K)
111REAL, DIMENSION(ngrid),        INTENT(IN)  :: wind1     ! wind speed first level (m/s)
112REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: pinterf ! pressure at interfaces(Pa)
113REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: x_in    ! variable at the beginning of timestep
114REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: K_in    ! eddy diffusivity coef at the beginning of time step
115REAL, DIMENSION(ngrid),        INTENT(IN)  :: x_surf  ! surface variable
116REAL, DIMENSION(ngrid),        INTENT(IN)  :: cdrag  ! drag coefficient
117REAL, DIMENSION(ngrid,nlay),   INTENT(OUT) :: x_predict  ! variable at the end of time step after explicit prediction
118
119
120integer i,k
121real ml,F1,rho
122real, dimension(ngrid) :: play1,temp1
123real, dimension(ngrid,nlay+1) :: K_big
124
125! computation of K_big
126
127play1(:)=play(:,1)
128temp1(:)=temp(:,1)
129
130! "big K" calculation
131do  k=2,nlay-1
132   do i=1,ngrid
133      rho=pinterf(i,k)/rd/(0.5*(temp(i,k-1)+temp(i,k)))
134      K_big(i,k)=rg*K_in(i,k)/(play(i,k)-play(i,k+1))*(rho**2)
135   enddo
136enddo
137! speficic treatment for k=nlay
138do i=1,ngrid
139   rho=pinterf(i,nlay)/rd/temp(i,nlay)
140   K_big(i,nlay)=rg*K_in(i,nlay)/(2*(play(i,nlay)-pinterf(i,nlay+1)))*(rho**2)
141enddo
142
143
144
145! x_predict calculation for 2<=k<=nlay-1
146do  k=2,nlay-1
147   do i=1,ngrid
148      ml=(pinterf(i,k)-pinterf(i,k+1))/rg
149      x_predict(i,k)=x_in(i,k)-dtime/ml*(-K_big(i,k+1)*x_in(i,k+1) &
150                  + (K_big(i,k)+K_big(i,k+1))*x_in(i,k) &
151                  - K_big(i,k)*x_in(i,k-1))
152   enddo
153enddo
154
155! Specific treatment for k=1
156do i=1,ngrid
157   ml=(pinterf(i,1)-pinterf(i,2))/rg
158   F1=-play1(i)/rd/temp1(i)*wind1(i)*cdrag(i)*(x_in(i,1)-x_surf(i)) ! attention convention sens du flux
159   x_predict(i,1)=x_in(i,1)-dtime/ml*(-K_big(i,2)*(x_in(i,2) - x_in(i,1))-F1)
160enddo
161
162! Specific treatment for k=nlay
163! flux at the top of the atmosphere=0
164do i=1,ngrid
165   ml=0.5*(pinterf(i,nlay)-pinterf(i,nlay+1))/rg
166   x_predict(i,nlay)=x_in(i,nlay)+dtime/ml*(K_big(i,nlay)*(x_in(i,nlay)-x_in(i,nlay-1)))
167enddo
168
169
170end subroutine atke_explicit_prediction
171
172
173
174end module lmdz_call_atke
Note: See TracBrowser for help on using the repository browser.