source: LMDZ5/trunk/libf/phymar/PHY_vdfTKE_RUN.f90 @ 3948

Last change on this file since 3948 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR


Integration of MAR physics

File size: 10.2 KB
RevLine 
[2089]1      subroutine PHY_vdfTKE_RUN
2
3!------------------------------------------------------------------------------+
4!                                                         Mon 17-Jun-2013  MAR |
5!     SubRoutine PHY_vdfTKE_RUN  includes Vertical Diffusion of TKE            |
6!                                                           and epsilon        |
7!                                                                              |
8!                                                                              |
9!     version 3.p.4.1 created by H. Gallee,               Tue 19-Mar-2013      |
10!           Last Modification by H. Gallee,               Mon 17-Jun-2013      |
11!                                                                              |
12!------------------------------------------------------------------------------+
13!                                                                              |
14!    INPUT: Kzm_AT           Vertical Turbulent Coeffic.(momentum) [m2/s2]     |
15!    ^^^^^^                                                                    |
16!                                                                              |
17!    INPUT / OUTPUT: The Vertical Turbulent Fluxes are included for:           |
18!    ^^^^^^^^^^^^^^                                                            |
19!         a) Turbulent Kinetic Energy             TKE_AT(_xyz)     [m2/s2]     |
20!         b) Turbulent Kinetic Energy Dissipation eps_AT(_xyz)     [m2/s3]     |
21!                                                                              |
22!   #OPTIONS: #De: Dirichlet Type Top Boundary Condit. for TKE_AT (TKE    )    |
23!   #^^^^^^^^                                            & eps_AT (epsilon)    |
24!                                                                              |
25!------------------------------------------------------------------------------+
26
27      use Mod_Real
28      use Mod_PHY____dat
29      use Mod_PHY____grd
30      use Mod_PHY_AT_grd
31      use Mod_PHY_AT_kkl
32      use Mod_PHY_DY_kkl
33
34
35
36! Local  Variables
37! ================
38
39      use Mod_vdfTKE_RUN
40
41
42      IMPLICIT NONE
43
44
45      real(kind=real8)                             ::  S3DSBC
46      real(kind=real8)                             ::  sige2k
47      real(kind=real8)                             ::  psa_sq
48! #De real(kind=real8)                             ::  TKEtop = 0.
49
50      integer                                      ::  i     ,j     ,k
51      integer                                      ::  k1           ,ikl
52
53
54
55
56!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
57!                                                         !
58! ALLOCATION
59! ==========
60
61      IF (it_RUN.EQ.1 .OR. FlagDALLOC)               THEN !
62          allocate ( S3D__1(mzp-1) )
63          allocate ( S3D__2(mzp-1) )
64          allocate ( S3D__3(mzp-1) )
65          allocate ( S3D__4(mzp-1) )
66          allocate ( S3D__5(mzp-1) )
67          allocate ( S3D__6(mzp-1) )
68          allocate ( S3D__7(mzp-1) )
69      END IF                                              !
70!                                                         !
71!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
72
73
74
75
76! INITIALIZATION
77! ==============
78
79      sige2k = sige / sigk
80
81
82
83!     =====================
84      DO      ikl = 1,kcolp
85!     =====================
86
87        i = ii__AP(ikl)
88        j = jj__AP(ikl)
89
90        psa_sq          = psa_DY(ikl)         *psa_DY(ikl)   
91
92
93
94
95! Vertical Diffusion of Turbulent Kinetic Energy
96! ==============================================
97
98
99! Tridiagonal Matrix Coefficients - TKE_AT
100! ----------------------------------------
101
102           k=mzp-1
103          S3DSBC        = -GravF2*0.5*(Kzm_AT(ikl,k)+Kzm_AT(ikl,k+1))  &!  SBC: TKE and epsilon, Atm SBL
104     &             * sigk *betaAT    * roamDY(ikl,k)*roa_DY(ikl,k  )   &!
105     &                   /(psa_sq    * dsigmi(k)    *dsigma(    k+1))
106          S3D__1(mzp-1) =  S3DSBC                                       !  SBC: TKE and epsilon, Atm SBL
107
108        DO k=mzp-2,1,-1
109          S3D__1(k)      =-GravF2*0.5*(Kzm_AT(ikl,k)+Kzm_AT(ikl,k+1))  &
110     &             * sigk *betaAT    * roamDY(ikl,k)*roa_DY(ikl,k  )   &
111     &                   /(psa_sq    * dsigmi(k)    *dsigma(    k+1))
112
113          S3D__3(k+1)    = S3D__1(k) * dsigmi(    k)/dsigmi(    k+1)   &
114     &                               / roa_DY(ikl,k)*roa_DY(ikl,k+1)
115        END DO
116
117          S3D__3(1)      = 0.0                                          !  UBC: Von Neuman     , Atm Top
118        DO k=   1,mzp-1
119          S3D__1(k) =      S3D__1(k) * dt__AT
120          S3D__3(k) =      S3D__3(k) * dt__AT
121          S3D__2(k) = 1.0 -S3D__3(k) -S3D__1(k)
122        END DO
123
124
125! Second Member of the Tridiagonal System - TKE_AT
126! ------------------------------------------------
127
128          S3D__4(1) =                                                   &
129     &    S3D__1(1) *a_b_AT*(TKE_AT(ikl,1)-TKE_AT(ikl,k1p(1)))
130! #De     S3D__1(1) = 0.0
131! #De     S3D__2(1) = 1.0
132! #De     S3D__4(1) = TKEtop
133
134        DO k=k1p(1),mzp-2
135          S3D__4(k) =                                                   &
136     &    S3D__1(k) *a_b_AT*(TKE_AT(ikl,    k )-TKE_AT(ikl,k1p(1)))     &
137     &   -S3D__3(k) *a_b_AT*(TKE_AT(ikl,k1m(k))-TKE_AT(ikl,    k ))
138        END DO
139
140           k=       mzp-1
141          S3D__4(k) =                                                   &
142     &    S3D__1(k)*(a_b_AT* TKE_AT(ikl,    k )-TKE_AT(ikl,k1p(1))      &
143     &                                         /betaAT            )     &
144     &   -S3D__3(k) *a_b_AT*(TKE_AT(ikl,k1m(k))-TKE_AT(ikl,    k ))
145
146          S3D__1(k) = 0.000
147
148!         S3D__4(mzp-1)=-(alphAT* TKE_AT(ikl,mzp-1)-TKE_AT(ikl,mzp  ))  &
149!    &          * GravF2 *0.5000*(Kzm_AT(ikl,mzp-1)+Kzm_AT(ikl,mzp-2))  &
150!    &                          * roamDY(ikl,mzp-1)*roamDY(ikl,mzp-1)   &
151!    &                 / (psa_sq* dsigmi(    mzp-1)*dsigma(    mzp-1))  &
152!    &   -S3D__3(mzp-1)*  a_b_AT*(TKE_AT(ikl,mzp-2)-TKE_AT(ikl,mzp-1))
153
154
155! Tridiagonal Matrix Inversion - TKE_AT
156! -------------------------------------
157
158             k1= 1
159! #De        k1= 2
160        DO k=k1,mzp-1
161          S3D__4(k) = S3D__4(k) +  TKE_AT(ikl,k)
162        END DO
163
164! Forward  Sweep
165! ~~~~~~~~~~~~~~
166          S3D__5(1) = S3D__2(1)
167          S3D__6(1) =-S3D__1(1)   /S3D__5(1)
168        DO   k=k1p(1),mzp-1
169          S3D__5(k) = S3D__3(k)   *S3D__6(k-1)+S3D__2(k)
170          S3D__6(k) =-S3D__1(k)   /S3D__5(k)
171        END DO
172          S3D__7(1) = S3D__4(1)   /S3D__5(1)
173        DO   k=k1p(1),mzp-1
174          S3D__7(k) =(S3D__4(k)   -S3D__3(k)                           &
175     &               *S3D__7(k-1))/S3D__5(k)
176        END DO
177
178! Backward Sweep
179! ~~~~~~~~~~~~~~
180        DO k=k1m(mzp-1),1,-1
181          S3D__7(k) = S3D__6(k)   *S3D__7(k+1)+S3D__7(k)
182        END DO
183
184
185        DO k=1,mzp-1
186          TrT_AT(ikl,k) =  TrT_AT(ikl,k)                               &
187     &              +(S3D__7(k)   -TKE_AT(ikl,k))  /dt__AT
188          TKE_AT(ikl,k) =                      S3D__7(k)
189        END DO
190
191
192! Vertical Diffusion of Dissipation
193! =================================
194
195
196! Update Tridiagonal Matrix Coefficients - eps_AT
197! -----------------------------------------------
198
199          S3D__1(mzp-1) =   S3DSBC    * dt__AT                          !  SBC: TKE and epsilon, Atm SBL
200        DO k=1,mzp-1
201          S3D__1(k) =       S3D__1(k) * sige2k
202          S3D__3(k) =       S3D__3(k) * sige2k
203          S3D__2(k) = 1.0 - S3D__3(k) - S3D__1(k)
204        END DO
205
206
207! Second Member of the Tridiagonal System - eps_AT
208! ------------------------------------------------
209
210          S3D__4(1) =                                                   &
211     &    S3D__1(1) *a_b_AT*(eps_AT(ikl,1)-eps_AT(ikl,k1p(1)))
212! #De     S3D__1(1) = 0.0
213! #De     S3D__2(1) = 1.0
214! #De     S3D__4(1) = eps_DI(i,j)
215
216        DO k=k1p(1),mzp-2
217          S3D__4(k) =                                                   &
218     &    S3D__1(k) *a_b_AT*(eps_AT(ikl,k)-eps_AT(ikl,k1p(k)))          &
219     &   -S3D__3(k) *a_b_AT*(eps_AT(ikl,k1m(k))-eps_AT(ikl,k))
220        END DO
221
222           k=       mzp-1
223          S3D__4(k) =                                                   &
224     &    S3D__1(k)*(a_b_AT* eps_AT(ikl,    k )-eps_AT(ikl,k1p(1))      &
225     &                                         /betaAT            )     &
226     &   -S3D__3(k) *a_b_AT*(eps_AT(ikl,k1m(k))-eps_AT(ikl,    k ))
227
228          S3D__1(k) = 0.000
229
230!         S3D__4(mzp-1)=-(alphAT* eps_AT(ikl,mzp-1)-eps_AT(ikl,mzp))    &
231!    &          * GravF2* 0.5000*(Kzm_AT(ikl,mzp-1)+Kzm_AT(ikl,mzp-2))  &
232!    &                          * roamDY(ikl,mzp-1)*roamDY(ikl,mzp-1)   &
233!    &                 / (psa_sq* dsigmi(    mzp-1)*dsigma(    mzp-1))  &
234!    &   -S3D__3(mzp-1)*  a_b_AT*(eps_AT(ikl,mzp-2)-eps_AT(ikl,mzp-1))
235
236
237! Tridiagonal Matrix Inversion - eps_AT
238! -------------------------------------
239
240             k1= 1
241! #De        k1= 2
242        DO k=k1,mzp-1
243          S3D__4(k)    = S3D__4(k) + eps_AT(ikl,k)
244        END DO
245
246! Forward  Sweep
247! ~~~~~~~~~~~~~~
248          S3D__5(1) = S3D__2(1)
249          S3D__6(1) =-S3D__1(1)   /S3D__5(1)
250        DO   k=k1p(1),mzp-1
251          S3D__5(k) = S3D__3(k)   *S3D__6(k-1)+S3D__2(k)
252          S3D__6(k) =-S3D__1(k)   /S3D__5(k)
253        END DO
254          S3D__7(1) = S3D__4(1)   /S3D__5(1)
255        DO   k=k1p(1),mzp-1
256          S3D__7(k) =(S3D__4(k)   -S3D__3(k)                            &
257     &               *S3D__7(k-1))/S3D__5(k)
258        END DO
259
260! Backward Sweep
261! ~~~~~~~~~~~~~~
262        DO k=k1m(mzp-1),1,-1
263          S3D__7(k) = S3D__6(k)   *S3D__7(k+1)+S3D__7(k)
264        END DO
265
266        DO k=1,mzp-1
267          eps_AT(ikl,k) =                      S3D__7(k)
268        END DO
269
270
271
272
273!     =====================
274      ENDDO ! ikl = 1,kcolp
275!     =====================
276
277
278
279
280!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281!                                                         !
282! DE-ALLOCATION
283! =============
284
285      IF (FlagDALLOC)                                THEN !
286          deallocate ( S3D__1 )
287          deallocate ( S3D__2 )
288          deallocate ( S3D__3 )
289          deallocate ( S3D__4 )
290          deallocate ( S3D__5 )
291          deallocate ( S3D__6 )
292          deallocate ( S3D__7 )
293      ENDIF
294!                                                         !
295!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
296
297
298
299
300      return
301      end subroutine PHY_vdfTKE_RUN
Note: See TracBrowser for help on using the repository browser.