source: LMDZ6/branches/contrails/libf/phylmd/ecrad/driver/test_spartacus_math.F90 @ 5418

Last change on this file since 5418 was 4773, checked in by idelkadi, 12 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (​https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


File size: 2.9 KB
Line 
1! (C) Copyright 2014- ECMWF.
2!
3! This software is licensed under the terms of the Apache Licence Version 2.0
4! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5!
6! In applying this licence, ECMWF does not waive the privileges and immunities
7! granted to it by virtue of its status as an intergovernmental organisation
8! nor does it submit to any jurisdiction.
9
10! This simple program tests the functioning of the matrix operations
11! in the radiation_matrix module
12
13program test_spartacus_math
14
15  use parkind1
16  use radiation_matrix
17
18  implicit none
19
20  integer, parameter :: n = 10, m = 3
21
22  real(jprb), dimension(n,m,m) :: A, B, C
23  real(jprb), dimension(n,m) :: v, w
24
25  ! Terms in exchange matrix
26  real(jprb), dimension(n) :: aa, bb, cc, dd
27
28  integer :: j
29
30  A = -1.0_jprb
31  A(:,1,1) = 10.0_jprb
32  A(:,2,2) = 5.0_jprb
33  A(:,2,1) = -0.5_jprb
34
35  B = 2.0_jprb
36  B(:,1,1) = -1.4_jprb
37  B(:,2,1) = 1.0e4_jprb
38
39  v(:,1) = 2.0_jprb
40  v(:,2) = 3.0_jprb
41  if (m == 3) then
42     v(:,m) = -1.5_jprb
43     A(:,m,m) = 4.5
44  end if
45
46  write(*,*) 'A ='
47  do j = 1,m
48     write(*,*) A(1,j,:)
49  end do
50
51  write(*,*) 'B ='
52  do j = 1,m
53     write(*,*) B(1,j,:)
54  end do
55
56  write(*,*) 'v = ', v(1,:)
57
58  C = mat_x_mat(n,n,m,A,B)
59  w = mat_x_vec(n,n,m,A,v)
60
61  write(*,*) 'C = A*B ='
62  do j = 1,m
63     write(*,*) C(1,j,:)
64  end do
65  write(*,*) 'w = A*v = ', w(1,:)
66
67  w = solve_vec(n,n,m,A,v)
68  C = solve_mat(n,n,m,A,B)
69
70
71  write(*,*) 'C = A\B ='
72  do j = 1,m
73     write(*,*) C(1,j,:)
74  end do
75  write(*,*) 'w = A\v = ', w(1,:)
76
77  call expm(n,n,m,A,IMatrixPatternDense)
78
79  write(*,*) 'expm(A) ='
80  do j = 1,m
81     write(*,*) A(1,j,:)
82  end do
83
84  ! Test fast_expm_exchange_3
85  aa = 2.0_jprb
86  bb = 3.0_jprb
87  cc = 5.0_jprb
88  dd = 7.0_jprb
89
90  if (m == 2) then
91
92    call fast_expm_exchange(n,n,aa,bb,A)
93   
94    write(*,*) 'fast_expm(A) = '
95    do j = 1,m
96      write(*,*) A(1,j,:)
97    end do
98
99    A = 0.0_jprb
100    A(:,1,1) = -aa
101    A(:,2,1) = aa
102    A(:,1,2) = bb
103    A(:,2,2) = -bb
104   
105    call expm(n,n,m,A,IMatrixPatternDense)
106   
107    write(*,*) 'expm(A) = '
108    do j = 1,m
109      write(*,*) A(1,j,:)
110    end do
111
112    ! Test zeros lead to identity matrix
113    aa = 0.0_jprb
114    call fast_expm_exchange(n,n,aa,aa,A)
115   
116    write(*,*) 'expm(zeros) = '
117    do j = 1,m
118      write(*,*) A(1,j,:)
119    end do
120
121  else
122   
123    call fast_expm_exchange(n,n,aa,bb,cc,dd,A)
124   
125    write(*,*) 'fast_expm(A) = '
126    do j = 1,m
127      write(*,*) A(1,j,:)
128    end do
129
130    A = 0.0_jprb
131    A(:,1,1) = -aa
132    A(:,2,1) = aa
133    A(:,1,2) = bb
134    A(:,2,2) = -bb-cc
135    A(:,3,2) = cc
136    A(:,2,3) = dd
137    A(:,3,3) = -dd
138   
139    call expm(n,n,m,A,IMatrixPatternDense)
140   
141    write(*,*) 'expm(A) = '
142    do j = 1,m
143      write(*,*) A(1,j,:)
144    end do
145
146    ! Test zeros lead to identity matrix
147    aa = 0.0_jprb
148    call fast_expm_exchange(n,n,aa,aa,aa,aa,A)
149   
150    write(*,*) 'expm(zeros) = '
151    do j = 1,m
152      write(*,*) A(1,j,:)
153    end do
154
155  end if
156   
157end program test_spartacus_math
Note: See TracBrowser for help on using the repository browser.