source: branches/iLoveclim/SOURCES/Temperature-routines/temp_col.f90 @ 77

Last change on this file since 77 was 77, checked in by dumas, 8 years ago

Merge branche iLOVECLIM sur rev 76
File size: 6.5 KB
Line 
1!> \file temp_col.f90
2!! Subroutines for vertical temperature calculation
3!<
4
5!> \namespace temp_col
6!!  Subroutines for vertical temperature calculation
7!<
8
9!> SUBROUTINE: Temp_column
10!! Subroutine for vertical temperature calculation
11!! Used modules:
12!! - use Icetemp_declar
13!! - use Tridiagmod
14!>
15
16Subroutine  Temp_col(Nz,Nzm,Nn,Newtempcol,Ctij,Flotij,Hij,Tbmerij,Ibase,Tij,Tpmpij,Aa,Bb,Cc,Rr,Hh,Ncond,Dee,Cm, &
17                                Dzm,Phidij,Ghfij,Ifail1,Tbdotij,DTT,Dou,Tsij,Bmeltij)
18
19!  Use Icetemp_declar
20  Use Tridiagmod
21 
22 
23  Implicit None
24  integer,intent(in) :: Nz,Nzm,Nn
25  Real,Dimension(Nz+Nzm),Intent(Out)  :: Newtempcol        !< Tableau 1d Vert. Des Temperatures En Sortie
26  Real,Dimension(Nz),Intent(In)       :: Ctij                !< Tableau 1d Vert. Conductivite Thermique
27  Logical,Intent(In)                  :: Flotij              !< Vrai Si Flotijtant (Test D'Archimede) 'O'
28  Real,Intent(In)                     :: Hij                 !< Ice Thickness  'O'
29  Real,Intent(In)                     :: Tbmerij             !< Temperature De La Mer A La Base De L'Ice Shelf
30  integer,intent(inout)               :: Ibase
31  real,dimension(Nz+Nzm),intent(inout) :: Tij             !< Temperature sur la colonne
32  real,dimension(Nz),intent(in)       :: Tpmpij               !< pressure melting point temperature in ice sheet 'o'
33  Real,Dimension(Nn),intent(inout) :: Aa    !< Work Arrays For Tridiag  !Dim Nn
34  Real,Dimension(Nn),intent(inout) :: Bb    !< Work Arrays For Tridiag  !Dim Nn
35  Real,Dimension(Nn),intent(inout) :: Cc    !< Work Arrays For Tridiag  !Dim Nn
36  Real,Dimension(Nn),intent(inout) :: Rr    !< Work Arrays For Tridiag  !Dim Nn
37  Real,Dimension(Nn),intent(inout) :: Hh    !< Work Arrays For Tridiag  !Dim Nn
38  Integer,intent(in)               :: Ncond !< Switch : Ncond=1 Thermal Conductivity In Mantle
39  real, intent(in) :: Dee
40  real,intent(in) :: Cm
41  real,intent(in) :: Dzm
42  real,intent(inout) :: Phidij     !< flux de chaleur lie a la deformation et glissement basal
43  real,intent(in) :: Ghfij      !< flux geothermique
44  integer,intent(out) :: Ifail1 !< Diagnostique erreur dans Tridiag
45  real, intent(out) :: Tbdotij  !< variation in time of basal temperature
46  real, intent(in) :: DTT
47  Real, intent(inout) :: Dou
48  real, intent(in) :: Tsij        !< surface ice temperature  'o'
49  real, intent(out) :: Bmeltij    !< fonte basale
50 
51 
52 
53  ! variables locales
54  integer :: K
55  Real :: Dzi,Tss
56  Real,Dimension(Nzm+1) :: Abis,Bbis,Cbis,Rbis,Hbis   
57
58  If (Hij.Gt.10.) Then
59
60     ! Conditions Aux Limites a La Base De La Glace
61     !________________________________________________
62
63     ! Base Froide
64     If ( ((Ibase.Eq.1).Or.(Ibase.Eq.4) &
65          .Or. ((Ibase.Eq.5).And.(Tij(Nz).Lt.Tpmpij(Nz)))) &
66          .And..Not.Flotij) Then
67
68        Ibase=1
69        Bb(Nz)=1.
70
71        If (Ncond.Eq.1) Then ! Avec Socle
72           Dzi=Hij*Dee*Cm/Ctij(Nz)
73           Aa(Nz)=-Dzm/(Dzm+Dzi)
74           Cc(Nz)=-Dzi/(Dzm+Dzi)
75           Rr(Nz)=Hij*Dee*Phidij*Dzm/Ctij(Nz)/(Dzm+Dzi)
76        Else                 ! Sans Socle
77           Aa(Nz)=-1.
78           Cc(Nz)=0.
79           Rr(Nz)=-(Ghfij-Phidij)/Ctij(Nz)*Hij*Dee
80        Endif
81
82     Else
83
84        ! Base Temperee Ou Shelf
85        ! ----------------------               
86        If (Ibase.Eq.5)  Ibase=2
87        Ibase=Max(Ibase,2)
88        Aa(Nz)=0.
89        Bb(Nz)=1.
90        Cc(Nz)=0.
91
92        If (.Not.Flotij) Then
93           Rr(Nz)=Tpmpij(Nz)
94        Else
95           Rr(Nz)=Tbmerij
96        End If
97
98     Endif
99
100
101     If (Ncond.Eq.1) Then ! Avec Socle
102        Do K=Nz+1,nz+nzm-1
103           Rr(K)=Tij(K)
104        End Do
105        Call Tridiag (Aa,Bb,Cc,Rr,Hh,nz+nzm,Ifail1)
106
107     Else                 ! Sans Socle
108
109        Call Tridiag (Aa,Bb,Cc,Rr,Hh,Nz,Ifail1)
110     Endif
111
112
113     ! Temperature Dans Le Socle Si Ncond=0, Lineaire Avec Le Gradient Ghf
114     If (Ncond.Eq.0) Then
115        Do K=Nz+1,nz+nzm
116           Hh(K)=Hh(Nz)-Dzm*(K-Nz)*Ghfij/Cm
117        End Do
118     Endif
119
120     Do K=1,nz+nzm
121        Newtempcol(K)=Hh(K)
122     End Do
123     Tbdotij=(Newtempcol(Nz)-Tij(Nz))/Dtt
124
125     !      ------------------------------------------------- Glace Tres Fine (H<10m Ou H=0)
126  Else If (Hij.Le.10.) Then
127
128     !       Pour Eviter Des Problemes Lorsque H Est Inferieur A 10 M
129     !       Profil Lineaire, Pente=Flux Geothermique
130     !       Le Cas Sans Glace Est Traite De La Meme Facon  (Dou=0)
131
132     ! ........................................ Avec Calcul Dans Le Socle
133     If (Ncond.Eq.1.And..Not.Flotij) Then   
134        If (Hij.Gt.0.) Then
135           ! Gradient Dans Le Socle
136           Dou=(Tij(Nz+1)-Tij(Nz))/Dzm*Cm     
137           Dou=Dou/Ctij(Nz)*Dee*Hij
138
139           Tss=Min(0.,Tsij)
140           Do K=1,Nz
141              Newtempcol(K)=Tss+Dou*(K-1.)
142           End Do
143
144        Else
145           ! Pas De Glace
146           Tss=Tsij
147           Do K=1,Nz
148              Newtempcol(K)=Tss
149           End Do
150        End If
151
152        ! Calcul Dans Le Socle Meme S'Il N'Y A Pas De Glace
153        Rr(Nz)=Newtempcol(Nz)
154        Aa(Nz)=0
155        Bb(Nz)=1
156        Cc(Nz)=0
157
158        Do K=Nz+1,nz+nzm-1
159           Rr(K)=Tij(K)
160        End Do
161
162        ! Creation De Nouveaux Tableaux Juste Pour Tridiag     
163        Do K=1,Nzm+1
164           Abis(K)=Aa(Nz-1+K)
165           Bbis(K)=Bb(Nz-1+K)
166           Cbis(K)=Cc(Nz-1+K)
167           Rbis(K)=Rr(Nz-1+K)
168        End Do
169
170        Call Tridiag (Abis,Bbis,Cbis,Rbis,Hbis,Nzm+1,Ifail1)
171
172        Do K=1,Nzm+1
173           Hh(Nz-1+K)=Hbis(K)
174        End Do
175
176        ! ........................................ Sans Calcul Dans Le Socle
177     Else
178        ! Calotte Posee                   
179        If ((Hij.Gt.0.).And..Not.Flotij) Then
180           Dou=-Ghfij/Ctij(Nz)*Dee*Hij
181           Tss=Min(0.,Tsij)
182           Do K=1,Nz
183              Newtempcol(K)=Tss+Dou*(K-1.)
184           End Do
185
186           !  Shelf
187        Else If ((Hij.Gt.0.).And.Flotij) Then
188           Tss=Min(0.,Tsij)
189           Dou=(Tbmerij-Tss)*Dee
190           Do K=1,Nz
191              Newtempcol(K)=Tss+Dou*(K-1.)
192           End Do
193
194           ! Pas De Glace
195        Else
196           Tss=Tsij     
197           Do K=1,Nz
198              Newtempcol(K)=Tss
199           End Do
200        End If
201
202        ! Temperature Dans Le Socle, Lineaire Avec Le Gradient Ghf
203        If (.Not.Flotij) Then
204           Do K=Nz+1,nz+nzm
205              Hh(K)=Newtempcol(Nz)-Dzm*(K-Nz)*Ghfij/Cm
206           End Do
207        Else
208           Do K=Nz+1,nz+nzm
209              Hh(K)=Tbmerij-Dzm*(K-Nz)*Ghfij/Cm
210           End Do
211        End If
212
213     Endif
214     ! Fin Du Test Socle
215
216
217     Do K=Nz+1,nz+nzm
218        Newtempcol(K)=Hh(K)
219     End Do
220
221     Tbdotij=(Newtempcol(Nz)-Tij(Nz))/Dtt
222     Bmeltij=0.
223     Ibase=5
224     Phidij=0.
225
226
227  End If
228
229End Subroutine Temp_col
Note: See TracBrowser for help on using the repository browser.