source: codes/icosagcm/trunk/src/timeloop_sw.f90 @ 12

Last change on this file since 12 was 12, checked in by ymipsl, 12 years ago

dynamico tree creation

YM
File size: 4.7 KB
Line 
1MODULE timeloop_sw_mod
2  USE genmod
3  USE transfert_mod
4  USE etat0_mod
5 
6  INTEGER,PARAMETER :: euler=1, leapfrog=2, leapfrog_matsuno=3, adam_bashforth=4
7 
8CONTAINS
9 
10  SUBROUTINE timeloop
11  USE field_mod
12  USE domain_mod
13!  USE wave_mod
14  USE caldyn_sw_mod
15  USE write_field
16  USE dimensions
17  USE geometry
18  USE transfert_mod
19  USE metric
20  USE dissip_mod
21  USE ioipsl
22
23  USE etat0_mod
24  IMPLICIT NONE
25  TYPE(t_field),POINTER :: f_h(:),f_hm1(:), f_hm2(:)
26  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
27  TYPE(t_field),POINTER :: f_dh(:),f_dhm1(:), f_dhm2(:)
28  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
29
30  REAL(rstd),POINTER :: h(:) ,hm1(:), hm2(:)
31  REAL(rstd),POINTER :: u(:) , um1(:), um2(:)
32  REAL(rstd),POINTER :: dh(:), dhm1(:), dhm2(:)
33  REAL(rstd),POINTER :: du(:), dum1(:), dum2(:)
34  REAL(rstd) :: dt
35  INTEGER :: ind
36  INTEGER :: it,i,j,n
37  INTEGER :: scheme
38  INTEGER :: matsuno_period
39  INTEGER :: itaumax
40 
41
42  dt=90.
43  CALL getin('dt',dt)
44
45  itaumax=100
46  CALL getin('itaumax',itaumax)
47 
48  scheme=leapfrog_matsuno
49  CALL getin('scheme',scheme)
50 
51  matsuno_period=5
52  CALL getin('matsuno_period',matsuno_period)
53  IF (scheme==leapfrog) matsuno_period=itaumax+1
54 
55 
56  CALL allocate_field(f_h,field_t,type_real)
57  CALL allocate_field(f_hm1,field_t,type_real)
58  CALL allocate_field(f_hm2,field_t,type_real)
59  CALL allocate_field(f_dh,field_t,type_real)
60  CALL allocate_field(f_dhm1,field_t,type_real)
61  CALL allocate_field(f_dhm2,field_t,type_real)
62
63  CALL allocate_field(f_u,field_u,type_real)
64  CALL allocate_field(f_um1,field_u,type_real)
65  CALL allocate_field(f_um2,field_u,type_real)
66  CALL allocate_field(f_du,field_u,type_real)
67  CALL allocate_field(f_dum1,field_u,type_real)
68  CALL allocate_field(f_dum2,field_u,type_real)
69
70
71  CALL allocate_caldyn
72 
73
74  CALL etat0_williamson(f_h,f_u)
75 
76  DO it=0,itaumax
77    PRINT *,"SW: It No :",It
78
79    CALL caldyn(f_h, f_u, f_dh, f_du)
80
81    IF (scheme==Euler) THEN
82      CALL  euler_scheme
83    ELSE IF (scheme==leapfrog) THEN
84      CALL leapfrog_scheme
85    ELSE IF (scheme==leapfrog_matsuno) THEN
86      CALL  leapfrog_matsuno_scheme
87    ELSE IF (scheme==adam_bashforth) THEN
88      CALL adam_bashforth_scheme
89    ENDIF
90
91  ENDDO
92 
93  CONTAINS
94
95    SUBROUTINE Euler_scheme
96    IMPLICIT NONE
97      INTEGER :: ind
98     
99      DO ind=1,ndomain
100        h=f_h(ind) ; u=f_u(ind) ; 
101        dh=f_dh(ind) ; du=f_du(ind) ;
102        h(:)=h(:)+dt*dh(:)
103        u(:)=u(:)+dt*du(:)
104      ENDDO
105     
106    END SUBROUTINE Euler_scheme
107   
108    SUBROUTINE leapfrog_scheme
109    IMPLICIT NONE
110      INTEGER :: ind
111
112      DO ind=1,ndomain
113        h=f_h(ind)   ; u=f_u(ind)   ; 
114        hm1=f_hm1(ind) ; um1=f_um1(ind) ; 
115        hm2=f_hm2(ind) ; um2=f_um2(ind) ; 
116        dh=f_dh(ind) ; du=f_du(ind) ;
117
118        IF (it==0) THEN
119          hm2(:)=h(:) ; um2(:)=u(:) 
120
121          h(:)=h(:)+dt*dh(:)
122          u(:)=u(:)+dt*du(:)
123          hm1(:)=h(:) ;  um1(:)=u(:) 
124
125        ELSE
126       
127          h(:)=hm2(:)+2*dt*dh(:)
128          u(:)=um2(:)+2*dt*du(:)
129
130          hm2(:)=hm1(:) ; um2(:)=um1(:) 
131          hm1(:)=h(:) ; um1(:)=u(:) 
132        ENDIF
133         
134      ENDDO
135    END SUBROUTINE leapfrog_scheme 
136 
137    SUBROUTINE leapfrog_matsuno_scheme
138    IMPLICIT NONE
139      INTEGER :: ind
140
141      DO ind=1,ndomain
142        h=f_h(ind)   ; u=f_u(ind) 
143        hm1=f_hm1(ind) ; um1=f_um1(ind) 
144        hm2=f_hm2(ind) ; um2=f_um2(ind) 
145        dh=f_dh(ind) ; du=f_du(ind) ; 
146
147       
148        IF (MOD(it,matsuno_period+1)==0) THEN
149          hm1(:)=h(:) ; um1(:)=u(:) 
150          h(:)=hm1(:)+dt*dh(:)
151          u(:)=um1(:)+dt*du(:)
152
153        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
154
155          h(:)=hm1(:)+dt*dh(:)
156          u(:)=um1(:)+dt*du(:)
157
158          hm2(:)=hm1(:) ;  um2(:)=um1(:) 
159          hm1(:)=h(:) ;  um1(:)=u(:) 
160
161        ELSE
162
163          h(:)=hm2(:)+2*dt*dh(:)
164          u(:)=um2(:)+2*dt*du(:)
165
166          hm2(:)=hm1(:) ; um2(:)=um1(:) 
167          hm1(:)=h(:) ;  um1(:)=u(:) 
168
169        ENDIF
170     
171      ENDDO
172     
173    END SUBROUTINE leapfrog_matsuno_scheme 
174         
175    SUBROUTINE adam_bashforth_scheme
176    IMPLICIT NONE
177      INTEGER :: ind
178
179      DO ind=1,ndomain
180        h=f_h(ind)   ; u=f_u(ind)   
181        dh=f_dh(ind) ; du=f_du(ind) 
182        dhm1=f_dhm1(ind) ; dum1=f_dum1(ind) 
183        dhm2=f_dhm2(ind) ; dum2=f_dum2(ind) 
184
185        IF (it==0) THEN
186          dhm1(:)=dh(:) ; dum1(:)=du(:) 
187          dhm2(:)=dhm1(:) ; dum2(:)=dum1(:)
188        ENDIF
189             
190        h(:)=h(:)+dt*(23*dh(:)-16*dhm1(:)+5*dhm2(:))/12
191        u(:)=u(:)+dt*(23*du(:)-16*dum1(:)+5*dum2(:))/12
192
193        dhm2(:)=dhm1(:) ; dum2(:)=dum1(:) 
194        dhm1(:)=dh(:) ; dum1(:)=du(:) 
195
196      ENDDO
197     
198     
199    END SUBROUTINE adam_bashforth_scheme
200
201  END SUBROUTINE timeloop
202   
203   
204 
205END MODULE timeloop_sw_mod
Note: See TracBrowser for help on using the repository browser.