Newer
Older
RODIER Quentin
committed
!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
MODULE MODE_TRIDIAG
IMPLICIT NONE
CONTAINS
SUBROUTINE TRIDIAG(KKA,KKU,KKL,PVARM,PA,PTSTEP,PEXPL,PIMPL, &
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
PRHODJ,PSOURCE,PVARP )
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! #################################################
!
!
!!**** *TRIDIAG* - routine to solve a time implicit scheme
!!
!!
!! PURPOSE
!! -------
! The purpose of this routine is to give a field PVARP at t+1, by
! solving an implicit tridiagonal system obtained by the
! discretization of the vertical turbulent diffusion. It should be noted
! that the degree of implicitness can be varied (PIMPL parameter) and the
! sources of evolution other than the turbulent diffusion can be taken
! into account through the PSOURCE field. PVARP is localized at a mass
! point.
!
!!** METHOD
!! ------
!! First, the Right Hand Side of the implicit equation is computed.
!! It is build as follows:
!! ZY = PVARM + PTSTEP*PSOURCE + DIFF_EXPLI
!! where PVARM is the variable at t-dt, PSOURCE the supplementary sources of
!! PVAR ( and not PVAR * PRHODJ !!) and DIFF_EXPLI is the explicit part
!! of the vertical turbulent diffusion. This operator is spatially
!! discretized as the implicit one, thus:
!! DIFF_EXPLI(k) = - PEXPL / PRHODJ(k) *
!! ( PA(k+1) * (PVARM(k+1) - PVARM(k) )
!! -PA(k) * (PVARM(k) - PVARM(k-1)) )
!! For the first level, only the upper part is considered, the lower one
!! is replaced by the turbulent surface flux (taken into account in the
!! PSOURCE(ikb) term).
!! DIFF_EXPLI(ikb) = - PEXPL / PRHODJ(ikb) *
!! ( PA(ikb+1) * (PVARM(ikb+1) - PVARM(ikb)) )
!! For the last level, only the lower part is considered, the upper one
!! is replaced by the turbulent flux which is taken equal to 0
!! (taken into account in the PSOURCE(ike) term).
!!
!! DIFF_EXPLI(ike) = + PEXPL / PRHODJ(ike) *
!! ( PA(ike) * (PVARM(ike) - PVARM(ike-1)) )
!!
!! Then, the classical tridiagonal algorithm is used to invert the
!! implicit operator. Its matrix is given by:
!!
!! ( b(ikb) c(ikb) 0 0 0 0 0 0 )
!! ( 0 a(ikb+1) b(ikb+1) c(ikb+1) 0 ... 0 0 0 )
!! ( 0 0 a(ikb+2) b(ikb+2) c(ikb+2). 0 0 0 )
!! .......................................................................
!! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 )
!! .......................................................................
!! ( 0 0 0 0 0 ...a(ike-1) b(ike-1) c(ike-1))
!! ( 0 0 0 0 0 ... 0 a(ike) b(ike) )
!!
!! ikb and ike represent the first and the last inner mass levels of the
!! model. The coefficients are:
!!
!! a(k) = PIMPL * PA(k)/PRHODJ(k)
!! b(k) = 1 - PIMPL * PA(k)/PRHODJ(k) - PIMPL * PA(k+1)/PRHODJ(k)
!! c(k) = PIMPL * PA(k+1)/PRHODJ(k)
!!
!! for all k /= ikb or ike
!!
!! b(ikb) = 1 - PIMPL * PA(ikb+1)/PRHODJ(ikb)
!! c(ikb) = PIMPL * PA(ikb+1)/PRHODJ(ikb)
!! (discretization of the upper part of the implicit operator)
!! b(ike) = 1 - PIMPL * PA(ike)/PRHODJ(ike)
!! a(ike) = PIMPL * PA(ike)/PRHODJ(ike)
!! (discretization of the lower part of the implicit operator)
!! Finally, the marginal points are prescribed.
!!
!! All these computations are purely vertical and vectorizations are
!! easely achieved by processing all the verticals in parallel.
!!
!! EXTERNAL
!! --------
!!
!! NONE
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! MODD_PARAMETERS
!! JPVEXT_TURB: number of vertical external points
!!
!! REFERENCE
!! ---------
!! Book 1 of Meso-NH documentation (chapter Turbulence)
!! Press et al: Numerical recipes (1986) Cambridge Univ. Press
!!
!! AUTHOR
!! ------
!! Joan Cuxart * INM and Meteo-France *
!!
!! MODIFICATIONS
!! -------------
!! Original August 29, 1994
!! Modification : January 29, 1995 Algorithm written with two
!! local variables less
!! (Cuxart, Stein) August 21, 1995 Bug correction for PRHODJ
!! (Stein) November 16, 1995 new version
!! (Stein) February 28, 1995 no inversion in the explicit case
!! (Seity) February 2012 add possibility to run with reversed
!! vertical levels
!! ---------------------------------------------------------------------
!
!* 0. DECLARATIONS
!
USE MODD_PARAMETERS
!
IMPLICIT NONE
!
!
!* 0.1 declarations of arguments
!
INTEGER, INTENT(IN) :: KKA !near ground array index
INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM ! variable at t-1
REAL, DIMENSION(:,:,:), INTENT(IN) :: PA ! upper diag. elements
REAL, INTENT(IN) :: PTSTEP ! Double time step
REAL, INTENT(IN) :: PEXPL,PIMPL ! weights of the temporal scheme
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! (dry rho)*J
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSOURCE ! source term of PVAR
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PVARP ! variable at t+1
!
!* 0.2 declarations of local variables
!
REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3)) :: ZY ,ZGAM
! RHS of the equation, 3D work array
REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZBET
! 2D work array
INTEGER :: JK ! loop counter
INTEGER :: IKB,IKE ! inner vertical limits
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
!
! ---------------------------------------------------------------------------
!
!* 1. COMPUTE THE RIGHT HAND SIDE
! ---------------------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TRIDIAG',0,ZHOOK_HANDLE)
IKTB=1+JPVEXT_TURB
IKT=SIZE(PVARM,3)
IKTE=IKT-JPVEXT_TURB
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
!
!
ZY(:,:,IKB) = PVARM(:,:,IKB) + PTSTEP*PSOURCE(:,:,IKB) - &
PEXPL / PRHODJ(:,:,IKB) * PA(:,:,IKB+KKL) * (PVARM(:,:,IKB+KKL) - PVARM(:,:,IKB))
!
DO JK=IKTB+1,IKTE-1
ZY(:,:,JK)= PVARM(:,:,JK) + PTSTEP*PSOURCE(:,:,JK) - &
PEXPL / PRHODJ(:,:,JK) * &
( PVARM(:,:,JK-KKL)*PA(:,:,JK) &
-PVARM(:,:,JK)*(PA(:,:,JK)+PA(:,:,JK+KKL)) &
+PVARM(:,:,JK+KKL)*PA(:,:,JK+KKL) &
)
END DO
!
ZY(:,:,IKE)= PVARM(:,:,IKE) + PTSTEP*PSOURCE(:,:,IKE) + &
PEXPL / PRHODJ(:,:,IKE) * PA(:,:,IKE) * (PVARM(:,:,IKE)-PVARM(:,:,IKE-KKL))
!
!
!* 2. INVERSION OF THE TRIDIAGONAL SYSTEM
! -----------------------------------
!
IF ( PIMPL > 1.E-10 ) THEN
!
!
! going up
!
ZBET(:,:) = 1. - PIMPL * PA(:,:,IKB+KKL) / PRHODJ(:,:,IKB) ! bet = b(ikb)
PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)
!
DO JK = IKB+KKL,IKE-KKL,KKL
ZGAM(:,:,JK) = PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK-KKL) / ZBET(:,:)
! gam(k) = c(k-1) / bet
ZBET(:,:) = 1. - PIMPL * ( PA(:,:,JK) * (1. + ZGAM(:,:,JK)) &
+ PA(:,:,JK+KKL) &
) / PRHODJ(:,:,JK)
! bet = b(k) - a(k)* gam(k)
PVARP(:,:,JK)= ( ZY(:,:,JK) - PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK) &
* PVARP(:,:,JK-KKL) &
) / ZBET(:,:)
! res(k) = (y(k) -a(k)*res(k-1))/ bet
END DO
! special treatment for the last level
ZGAM(:,:,IKE) = PIMPL * PA(:,:,IKE) / PRHODJ(:,:,IKE-KKL) / ZBET(:,:)
! gam(k) = c(k-1) / bet
ZBET(:,:) = 1. - PIMPL * ( PA(:,:,IKE) * (1. + ZGAM(:,:,IKE)) &
) / PRHODJ(:,:,IKE)
! bet = b(k) - a(k)* gam(k)
PVARP(:,:,IKE)= ( ZY(:,:,IKE) - PIMPL * PA(:,:,IKE) / PRHODJ(:,:,IKE) &
* PVARP(:,:,IKE-KKL) &
) / ZBET(:,:)
! res(k) = (y(k) -a(k)*res(k-1))/ bet
!
! going down
!
DO JK = IKE-KKL,IKB,-1*KKL
PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+KKL) * PVARP(:,:,JK+KKL)
END DO
!
ELSE
!
PVARP(:,:,IKTB:IKTE) = ZY(:,:,IKTB:IKTE)
!
END IF
!
!
!* 3. FILL THE UPPER AND LOWER EXTERNAL VALUES
! ----------------------------------------
!
PVARP(:,:,KKA)=PVARP(:,:,IKB)
PVARP(:,:,KKU)=PVARP(:,:,IKE)
!
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TRIDIAG',1,ZHOOK_HANDLE)
END SUBROUTINE TRIDIAG
END MODULE MODE_TRIDIAG