source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/call_atke_mod.F90 @ 5134

Last change on this file since 5134 was 4669, checked in by Laurent Fairhead, 16 months ago

Merged with trunk revision 4586 corresponding to june 2023 testing
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, &
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, &
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.