source: trunk/SOURCES/Temperature-routines/advec_icetemp.f90 @ 111

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

OpenMP parallelization in Temperature-routines, deformation_mod_2lois and resol_adv_diff_2D-sept2009
File size: 4.3 KB
Line 
1!> \file advec_icetemp.f90
2!! Subroutines for temperature advection in the heat transfer equation
3!<
4
5!> \namespace advec_icetemp
6!! Subroutines to calculate the temperature advection
7!<
8
9!> SUBROUTINE: Advec_icetemp
10!! Implementation of Routines to calculate advection in the heat transfer equation
11!! for a column
12!! Used modules:
13!!    - use Icetemp_declar
14!>
15
16Subroutine  Advec_icetemp(Nz,Nzm,Nn,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m,Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n, &
17        Uxij,Uxi_plus_un,Uyij,Uyj_plus_un,Tij,Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un,Hij,Tsij,Aa,Bb,Cc,Rr, &
18        Dx11,Dou,DTT,Dee,Uzrij,Dttdx)
19
20!  use icetemp_declar
21
22  Implicit None
23    !< Arguments
24  Integer,                     intent(in)  ::  Nz,Nzm,Nn      !< Taille des tableaux
25  Real, Dimension(Nz),         intent(in)   :: Prodq_m        !< Tableau 1d Vert. heat production
26  Real, Dimension(Nz),         intent(in)   :: Cro_m          !< Tableau 1d Vert. heat capcity
27  Real, Dimension(Nz),         intent(inout):: Advecx_m       !< Advection selon x
28  Real, Dimension(Nz),         intent(inout):: Advecy_m       !< Advection selon y
29  Real, Dimension(Nz),         intent(inout):: Advec_m        !< Advection total
30  Real, Dimension(Nz),         intent(in)   :: Ct_m           !< Tableau 1d Vert.  thermal cond.
31 
32  Integer, intent(in) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n
33  real, dimension(Nz), Intent(in) :: Uxij, Uxi_plus_un
34  real, dimension(Nz), Intent(in) :: Uyij, Uyj_plus_un
35  real, dimension(Nz+Nzm), Intent(inout) :: Tij
36  real, dimension(Nz+Nzm), Intent(in) :: Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un
37  real, intent(in) :: Hij
38  real,intent(in) :: Tsij
39  Real,Dimension(Nn),intent(inout) :: Aa    !< Work Arrays For Tridiag  !Dim Nn
40  Real,Dimension(Nn),intent(inout) :: Bb    !< Work Arrays For Tridiag  !Dim Nn
41  Real,Dimension(Nn),intent(inout) :: Cc    !< Work Arrays For Tridiag  !Dim Nn
42  Real,Dimension(Nn),intent(inout) :: Rr    !< Work Arrays For Tridiag  !Dim Nn
43  Real, intent(in) :: Dx11
44  Real, intent(inout) :: Dou
45  real, intent(in) :: DTT
46  real, intent(in) :: Dee
47  real, dimension(Nz) :: Uzrij
48  real :: Dttdx
49 
50
51 
52  ! variables locales
53  Integer :: K
54  Real :: Dzz,Dah,Ct_haut,Ct_bas
55
56 
57  Do K=2,Nz-1
58
59
60     ! On Remplit Deux Tableaux Pour eviter Des Problemes D'Asymetrie
61     !         Avection Selon X  (advecx_m en general > 0)
62     !         ----------------
63     Advecx_m(K) = Iadvec_w * Uxij(K) *          & ! ux west if upwind
64          (Tij(K) - Ti_moins_un(K))  & ! west T gradient
65           
66          + Iadvec_e * Uxi_plus_un(K)*      & ! ux east if upwind
67          (Ti_plus_un(K)-Tij(K))     ! east T gradient
68
69     Advecx_m(K) =  Advecx_m(K) * Dx11                                        !Dx11=1/Dx
70
71
72     !         Avection Selon Y
73     !         ----------------
74     Advecy_m(K) =   Iadvec_s * Uyij(K) *        & ! uy sud if upwind
75          (Tij(K) - Tj_moins_un(K))& ! south T gradient
76         
77          + Iadvec_n * Uyj_plus_un(K) *   & ! uy nord is upwind
78          (Tj_plus_un(K)-Tij(K))    ! north T gradient
79
80     Advecy_m(K) =  Advecy_m(K) * Dx11
81
82     Advec_m(K)=Advecx_m(K)+Advecy_m(K)  ! Advection Totale
83
84  End Do
85
86  !     Fin Boucle K, Ecriture Matrice Tridiag Dans La Glace
87
88
89  !      -----------------------------Cas General (H>10m)
90  thick_ice:  If (Hij.Gt.10.) Then
91
92     !        Variables De Calcul Dans La Glace
93     ! dou = dtt/dz^2
94     Dou=Dtt/Dee/Dee/ Hij / Hij 
95
96     ! Dah = dtt/dz
97     Dah=Dtt /Dee /Hij
98
99     ! thermal conductivity at mid point just below the surface                   
100     Ct_bas=2.*(Ct_m(1)*Ct_m(2))/(Ct_m(1)+Ct_m(2))           
101
102     ! surface temperature : cannot go above 0 celsius
103     Tij(1)=Min(0.,Tsij)
104
105     ! surface boundary condition
106     Aa(1)=0.
107     Bb(1)=1.
108     Cc(1)=0.
109     Rr(1)=Tij(1)
110
111     Do K=2,Nz-1
112        Dzz=Dou/Cro_m(K)
113
114        ! Conductivite Au Milieu Des Mailles
115        Ct_haut=Ct_bas
116        Ct_bas=2.*(Ct_m(K)*Ct_m(K+1))/(Ct_m(K)+Ct_m(K+1))
117
118        ! Advection Verticale Centree
119        Aa(K) = -Dzz * Ct_haut - Uzrij(K) * Dah / 2.  ! lower diag
120        Cc(K) = -Dzz * Ct_bas  + Uzrij(K) * Dah / 2.  ! upper diag
121        Bb(K) =  1.+ Dzz * (Ct_haut+Ct_bas)                                ! diag
122       
123        Rr(K) =  Tij(K) + Prodq_m(K) * Dtt   ! vector
124        Rr(K) =  Rr(K)- Dttdx * (Advec_m(K))
125
126     End DO
127  End If thick_ice
128
129End Subroutine Advec_icetemp
130
Note: See TracBrowser for help on using the repository browser.