SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
tridiag_surf.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! #########
6  SUBROUTINE tridiag_surf(PVARM,PF,PDFDDTDZ,PEXT,PDEXTDV,PTSTEP, &
7  pdzz,pdzm,pvarp,oimpl,palfa,pbeta )
8 ! #################################################
9 !
10 !
11 !!**** *TRIDIAG_SURF* - routine to solve a time implicit scheme
12 !!
13 !!
14 !! PURPOSE
15 !! -------
16 ! The purpose of this routine is to give a field PVARP at t+1, by
17 ! solving an implicit TRIDIAGonal system obtained by the
18 ! discretization of the vertical turbulent diffusion. It should be noted
19 ! that the degree of implicitness can be varied (ZIMPL parameter) and that
20 ! the function of F(dT/dz) must have been linearized.
21 ! PVARP is localized at a mass point.
22 !
23 !!** METHOD
24 !! ------
25 !!
26 !! [T(+) - T(-)]/Dt = -d{ F + dF/d(dT/dz) * impl * [dT/dz(+) - dT/dz(-)] }/dz
27 !! + Ext + dExt/dT * impl [ T(+) - T(-)]
28 !!
29 !! It is discretized as follows:
30 !!
31 !! PDZM(k)*PVARP(k)/PTSTEP
32 !! =
33 !! PDZM(k)*PVARM(k)/PTSTEP
34 !! - (PDZM(k+1)+PDZM(k) )/2. * PF(k+1)/PDZZ(k+1)
35 !! + (PDZM(k) +PDZM(k-1))/2. * PF(k) /PDZZ(k)
36 !! + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2
37 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARP(k+1)/PDZZ(k+1)**2
38 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2
39 !! + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARP(k) /PDZZ(k+1)**2
40 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2
41 !! + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARP(k) /PDZZ(k)**2
42 !! + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2
43 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARP(k-1)/PDZZ(k)**2
44 !! + PDZM(k) * PEXT(k)
45 !! + PDZM(k) * ZIMPL* PDEXTDV(k) * PVARP(k)
46 !! - PDZM(k) * ZIMPL* PDEXTDV(k) * PVARM(k)
47 !!
48 !!
49 !!
50 !! The system to solve is:
51 !!
52 !! A*PVARP(k-1) + B*PVARP(k) + C*PVARP(k+1) = Y(k)
53 !!
54 !!
55 !! The RHS of the linear system in PVARP writes:
56 !!
57 !! y(k) = PDZM(k)*PVARM(k)/PTSTEP
58 !! - (PDZM(k+1)+PDZM(k) )/2. * PF(k+1)/PDZZ(k+1)
59 !! + (PDZM(k) +PDZM(k-1))/2. * PF(k) /PDZZ(k)
60 !! + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2
61 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2
62 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2
63 !! + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2
64 !! + PDZM(k) * PEXT(k)
65 !! - PDZM(k) * PDEXTDV(k) * PVARM(k)
66 !!
67 !!
68 !! Then, the classical TRIDIAGonal algorithm is used to invert the
69 !! implicit operator. Its matrix is given by:
70 !!
71 !! ( b(ikb) c(ikb) 0 0 0 0 0 0 )
72 !! ( 0 a(ikb+1) b(ikb+1) c(ikb+1) 0 ... 0 0 0 )
73 !! ( 0 0 a(ikb+2) b(ikb+2) c(ikb+2). 0 0 0 )
74 !! .......................................................................
75 !! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 )
76 !! .......................................................................
77 !! ( 0 0 0 0 0 ...a(ike-1) b(ike-1) c(ike-1))
78 !! ( 0 0 0 0 0 ... 0 a(ike) b(ike) )
79 !!
80 !! ikb and ike represent the first and the last inner mass levels of the
81 !! model. The coefficients are:
82 !!
83 !! a(k) = + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) /PDZZ(k)**2
84 !! b(k) = PDZM(k) / PTSTEP
85 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2
86 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) /PDZZ(k)**2
87 !! - PDZM(k) * PDEXTDV(k)
88 !! c(k) = + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2
89 !!
90 !! for all k /= ikb or ike
91 !!
92 !!
93 !! b(ikb) = PDZM(ikb) / PTSTEP
94 !! -(PDZM(ikb+1)+PDZM(ikb))/2.*ZIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2
95 !! c(ikb) = +(PDZM(ikb+1)+PDZM(ikb))/2.*ZIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2
96 !!
97 !! b(ike) = PDZM(ike) / PTSTEP
98 !! -(PDZM(ike)+PDZM(ike-1))/2.*ZIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2
99 !! a(ike) = +(PDZM(ike)+PDZM(ike-1))/2.*ZIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2
100 !!
101 !!
102 !! EXTERNAL
103 !! --------
104 !!
105 !! NONE
106 !!
107 !! IMPLICIT ARGUMENTS
108 !! ------------------
109 !!
110 !! REFERENCE
111 !! ---------
112 !! Press et al: Numerical recipes (1986) Cambridge Univ. Press
113 !!
114 !! AUTHOR
115 !! ------
116 !! V. Masson * Meteo-France *
117 !!
118 !! MODIFICATIONS
119 !! -------------
120 !! Original 04/2003 (from tridiag.f90)
121 !! 09/2007 implicitation with surface flux at lowest level
122 !! ---------------------------------------------------------------------
123 !
124 !* 0. DECLARATIONS
125 !
126 !
127 USE yomhook ,ONLY : lhook, dr_hook
128 USE parkind1 ,ONLY : jprb
129 !
130 IMPLICIT NONE
131 !
132 !
133 !* 0.1 declarations of arguments
134 !
135 REAL, DIMENSION(:,:), INTENT(IN) :: pvarm ! variable at t-1 at mass point
136 REAL, DIMENSION(:,:), INTENT(IN) :: pf ! flux in dV/dt=-dF/dz at flux point
137 REAL, DIMENSION(:,:), INTENT(IN) :: pdfddtdz! dF/d(dV/dz) at flux point
138 REAL, DIMENSION(:,:), INTENT(IN) :: pext ! External forces in dT/dt=Ext
139 ! ! (Except surf. flux) at mass point
140 REAL, DIMENSION(:,:), INTENT(IN) :: pdextdv ! dExt/dV
141 ! ! at mass point
142 REAL, INTENT(IN) :: ptstep ! time step
143 REAL, DIMENSION(:,:), INTENT(IN) :: pdzz ! Dz at flux point
144 REAL, DIMENSION(:,:), INTENT(IN) :: pdzm ! Dz at mass point
145 !
146 REAL, DIMENSION(:,:), INTENT(OUT):: pvarp ! variable at t+1 at mass point
147 LOGICAL, OPTIONAL, INTENT(IN) :: oimpl ! true if implicit coupling
148 ! ! information between lvl 1 and
149 ! ! surface flux scheme must be
150 ! ! returned.
151 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: palfa ! V+(1) = alfa F(1) + beta
152 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pbeta ! V+(1) = alfa F(1) + beta
153 !
154 !
155 !* 0.2 declarations of local variables
156 !
157 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: zvarp
158 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: zdz_dfddtdz_o_dz2
159 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: za, zb, zc
160 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: zy ,zgam
161  ! RHS of the equation, 3D work array
162 REAL, DIMENSION(SIZE(PVARM,1)) :: zbet
163  ! 2D work array
164 INTEGER :: jk ! loop counter
165 INTEGER :: ik ! vertical limits
166 !
167 LOGICAL :: gimpl
168 REAL :: zimpl ! implicitation coefficient
169 REAL(KIND=JPRB) :: zhook_handle
170 ! ! for SBL scheme solving
171 ! ---------------------------------------------------------------------------
172 !
173 !* 1. Preliminaries
174 ! -------------
175 !
176 IF (lhook) CALL dr_hook('TRIDIAG_SURF',0,zhook_handle)
177 zimpl = 1
178 !
179 ik = SIZE(pf,2)
180 !
181 zdz_dfddtdz_o_dz2 = pdfddtdz/pdzz
182 !
183 za=0.
184 zb=0.
185 zc=0.
186 zy=0.
187 !
188 IF (present(oimpl)) THEN
189  gimpl = oimpl
190 ELSE
191  gimpl = .false.
192 END IF
193 !
194 !* 2. COMPUTE THE RIGHT HAND SIDE
195 ! ---------------------------
196 !
197 zy(:,1) = pdzm(:,1)*pvarm(:,1)/ptstep &
198  - pf(:,2) &
199  + pf(:,1) &
200  - zdz_dfddtdz_o_dz2(:,1) * zimpl * pvarm(:,1)&
201  + zdz_dfddtdz_o_dz2(:,2) * zimpl * pvarm(:,2)&
202  - zdz_dfddtdz_o_dz2(:,2) * zimpl * pvarm(:,1)&
203  + pdzm(:,1)*pext(:,1) &
204  - pdzm(:,1)*pdextdv(:,1) * zimpl * pvarm(:,1)
205 !
206 DO jk=2,ik-1
207  zy(:,jk) = pdzm(:,jk)*pvarm(:,jk)/ptstep &
208  - pf(:,jk+1) &
209  + pf(:,jk ) &
210  + zdz_dfddtdz_o_dz2(:,jk+1) * zimpl * pvarm(:,jk+1) &
211  - zdz_dfddtdz_o_dz2(:,jk+1) * zimpl * pvarm(:,jk ) &
212  - zdz_dfddtdz_o_dz2(:,jk ) * zimpl * pvarm(:,jk ) &
213  + zdz_dfddtdz_o_dz2(:,jk ) * zimpl * pvarm(:,jk-1) &
214  + pdzm(:,jk)*pext(:,jk) &
215  - pdzm(:,jk)*pdextdv(:,jk) * zimpl * pvarm(:,jk)
216 END DO
217 !
218 !* upper level is fixed (atmospheric forcing)
219 ! turbulent flux divergence is supposed equal to ZERO (inertial sublayer).
220 ! other terms are kept
221 !
222 zy(:,ik) = pdzm(:,ik)*pvarm(:,ik)/ptstep &
223  + pdzm(:,ik)*pext(:,ik) &
224  - pdzm(:,ik)*pdextdv(:,ik) * zimpl * pvarm(:,ik)
225 !
226 !
227 !* 3. INVERSION OF THE TRIDIAGONAL SYSTEM
228 ! -----------------------------------
229 !
230 !
231 !* 3.1 arrays A, B, C
232 ! --------------
233 !
234  zb(:,1) = pdzm(:,1)/ptstep &
235  - zdz_dfddtdz_o_dz2(:,1) * zimpl &
236  - zdz_dfddtdz_o_dz2(:,2) * zimpl &
237  - pdzm(:,1)*pdextdv(:,1)
238  zc(:,1) = zdz_dfddtdz_o_dz2(:,2) * zimpl
239 
240  DO jk=2,ik-1
241  za(:,jk) = zdz_dfddtdz_o_dz2(:,jk ) * zimpl
242  zb(:,jk) = pdzm(:,jk)/ptstep &
243  - zdz_dfddtdz_o_dz2(:,jk+1) * zimpl &
244  - zdz_dfddtdz_o_dz2(:,jk ) * zimpl &
245  - pdzm(:,jk)*pdextdv(:,jk)
246  zc(:,jk) = zdz_dfddtdz_o_dz2(:,jk+1) * zimpl
247  END DO
248 
249  !* upper level is fixed (atmospheric forcing)
250  ! turbulent flux divergence is supposed equal to ZERO (inertial sublayer).
251  ! other terms are kept
252  za(:,ik) = 0.
253  zb(:,ik) = pdzm(:,ik)/ptstep &
254  - pdzm(:,ik)*pdextdv(:,ik)
255 !
256 !* 3.2 going down
257 ! ----------
258 !
259  zbet(:) = zb(:,ik) ! bet = b(ik)
260  zvarp(:,ik) = zy(:,ik) / zbet(:)
261 
262  !
263  DO jk = ik-1,2,-1
264  zgam(:,jk) = za(:,jk+1) / zbet(:)
265  ! gam(k) = c(k-1) / bet
266  zbet(:) = zb(:,jk) - zc(:,jk) * zgam(:,jk)
267  ! bet = b(k) - a(k)* gam(k)
268  zvarp(:,jk)= ( zy(:,jk) - zc(:,jk) * zvarp(:,jk+1) ) / zbet(:)
269  ! res(k) = (y(k) -a(k)*res(k-1))/ bet
270  END DO
271  ! special treatment for the lowest level
272  zgam(:,1) = za(:,2) / zbet(:)
273  ! gam(k) = c(k-1) / bet
274  zbet(:) = zb(:,1) - zc(:,1) * zgam(:,1)
275  ! bet = b(k) - a(k)* gam(k)
276  zvarp(:,1)= ( zy(:,1) - zc(:,1) * zvarp(:,2) ) / zbet(:)
277  ! res(k) = (y(k) -a(k)*res(k-1))/ bet
278 !
279 !
280 !* 3.3 going up
281 ! --------
282 !
283  DO jk = 2,ik
284  zvarp(:,jk) = zvarp(:,jk) - zgam(:,jk-1) * zvarp(:,jk-1)
285  END DO
286 
287 !
288 !
289 !* 3.4 updates prognostic variable
290 ! ---------------------------
291 !
292 IF (.NOT. gimpl) THEN
293  pvarp=zvarp
294 ELSE
295  palfa=1./zbet(:)
296  pbeta=zvarp(:,1)
297  pvarp=pvarm
298 END IF
299 IF (lhook) CALL dr_hook('TRIDIAG_SURF',1,zhook_handle)
300 !
301 !-------------------------------------------------------------------------------
302 !
303 END SUBROUTINE tridiag_surf
tridiag_surf
subroutine tridiag_surf(PVARM, PF, PDFDDTDZ, PEXT, PDEXTDV, PTSTEP, PDZZ, PDZM, PVARP, OIMPL, PALFA, PBETA)
Definition: tridiag_surf.F90:6
Generated on Fri May 13 2016 09:31:05 for SURFEX by   doxygen 1.8.5