-
Notifications
You must be signed in to change notification settings - Fork 19
/
matinv.F90
163 lines (159 loc) · 5.27 KB
/
matinv.F90
1
2
3
4
5
6
7
8
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
! ------------------------------------------------------------------
! --- matrix inversion subroutines for implicit solution of vertical
! --- diffusion equation - tri-diagonal matrix
! ------------------------------------------------------------------
!
subroutine tridcof(diff,tri,nlayer,tcu,tcc,tcl)
use mod_xc ! HYCOM communication interface
use mod_cb_arrays ! HYCOM saved arrays
implicit none
!
! --- compute coefficients for tridiagonal matrix (dimension=kdm).
! --- Note: tcu(1) = 0. and tcl(kdm+1) = 0. are necessary conditions.
!
! --- input
real diff(kdm+1) ! diffusivity profile on interfaces
integer nlayer
!
! --- output
real tcu(kdm), & ! upper coeff. for (k-1) on k line of trid.matrix
tcc(kdm), & ! central ... (k ) ..
tcl(kdm) ! lower ..... (k-1) ..
! --- common tridiagonal factors
real tri(kdm,0:1) ! dt/dz/dz factors in trid. matrix
!
! --- local
integer k
!
! --- in the surface layer
tcu(1)=0.
tcc(1)=1.+tri(1,1)*diff(2) ! 1.+ delt1/h(1)/dzb(1)*diff(2)
tcl(1)= -tri(1,1)*diff(2) ! - delt1/h(1)/dzb(1)*diff(2)
!
! --- inside the domain
do 10 k=2,nlayer
tcu(k)= -tri(k,0)*diff(k )
tcc(k)=1.+tri(k,1)*diff(k+1)+tri(k,0)*diff(k)
tcl(k)= -tri(k,1)*diff(k+1)
10 continue
!
! --- in the bottom layer
tcl(nlayer)= 0.
return
end
!**********************************************************************
subroutine tridrhs(h,yo,diff,ghat,ghatflux,tri,nlayer,rhs)
use mod_xc ! HYCOM communication interface
use mod_cb_arrays ! HYCOM saved arrays
implicit none
!
! --- compute right hand side of tridiagonal matrix for scalar fields:
! --- = yo (old field)
! --- + flux-divergence of ghat
! --- + flux-divergence of non-turbulant fluxes
!
! --- note: if surface and bottom fluxes are nonzero, the following must apply
! --- sfc. lyr. needs +delt1/h(1)*surfaceflux
! --- bot. lyr. needs +delt1/h(nlayer)*diff(nlayer+1)/
! --- dzb(nlayer)*yo(nlayer+1)
!
! --- input
real h(kdm), & ! layer thickness
yo(kdm+1), & ! old profile
diff(kdm+1), & ! diffusivity profile on interfaces
ghat(kdm+1), & ! ghat turbulent flux
ghatflux ! surface flux for ghat: includes solar flux
integer nlayer
!
! --- output
real rhs(kdm) ! right hand side
!
real tri(kdm,0:1) ! dt/dz/dz factors in trid. matrix
!
! --- local
integer k
!
! --- in the top layer
rhs(1)=yo(1)+delt1/h(1)*(ghatflux*diff(2)*ghat(2))
!
! --- inside the domain
do 10 k=2,nlayer-1
rhs(k)=yo(k)+delt1/h(k)* &
(ghatflux*(diff(k+1)*ghat(k+1)-diff(k)*ghat(k)))
10 continue
!
! --- in the bottom layer
k=nlayer
rhs(k)=yo(k)+delt1/h(k)* &
(ghatflux*(diff(k+1)*ghat(k+1)-diff(k)*ghat(k)))
!
return
end
!**********************************************************************
subroutine tridmat(tcu,tcc,tcl,nlayer,h,rhs,yo,yn,diff, i,j)
use mod_xc ! HYCOM communication interface
use mod_cb_arrays ! HYCOM saved arrays
implicit none
!
! --- solve tridiagonal matrix for new vector yn, given right hand side
! --- vector rhs.
!
! --- note: if surface and bottom fluxes are nonzero, the following must apply
! --- surface layer needs +delt1*surfaceflux/(h(1)*bet)
! --- bottom layer needs +tri(nlayer,1)*diff(nlayer+1)*yo(nlayer+1))/bet
!
! --- input
real tcu (kdm), & ! upper coeff. for (k-1) on k line of tridmatrix
tcc (kdm), & ! central ... (k ) ..
tcl (kdm), & ! lower ..... (k-1) ..
h (kdm), & ! layer thickness
rhs(kdm), & ! right hand side
yo(kdm+1), & ! old field
diff(kdm+1), & ! diffusivity profile
gam(kdm) ! temporary array for tridiagonal solver
real bet ! ...
integer nlayer, & ! number of active layers <=kdm
i,j ! local grid point
!
! --- output
real yn(kdm+1) ! new field
!
! --- local
integer k
!
! --- solve tridiagonal matrix.
bet=tcc(1)
yn(1)=rhs(1)/bet ! surface
do 21 k=2,nlayer
gam(k)=tcl(k-1)/bet
bet=tcc(k)-tcu(k)*gam(k)
if(bet.eq.0.) then
write(lp,*)
write(lp,*) '** algorithm for solving tridiagonal matrix fails'
write(lp,*) '** i,j=',i0+i,j0+j !global grid point
write(lp,*) '** bet=',bet
write(lp,*) '** k=',k,' tcc=',tcc(k),' tcu=',tcu(k), &
' gam=',gam(k)
call flush(lp)
call xchalt('(tridmat)')
stop '(tridmat)'
! bet=1.E-12
endif
yn(k) = (rhs(k) - tcu(k) *yn(k-1) )/bet
! to avoid "Underflow" at single precision on the sun
! yni = (rhs(k) - tcu(k) *yn(k-1) )/bet
! if(yni.lt.0.) then
! yn(k) =min( (rhs(k) - tcu(k) *yn(k-1) )/bet ,-1.E-12 )
! else
! yn(k) = max( (rhs(k) - tcu(k) *yn(k-1) )/bet , 1.E-12 )
! endif
21 continue
!
do 22 k=nlayer-1,1,-1
yn(k)=yn(k)-gam(k+1)*yn(k+1)
22 continue
!
yn(nlayer+1)=yo(nlayer+1)
!
return
end