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

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

Update on 3D dynamic

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_sw_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  CALL init_dissip(dt)
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
82    IF (scheme==Euler) THEN
83      CALL  euler_scheme
84    ELSE IF (scheme==leapfrog) THEN
85      CALL leapfrog_scheme
86    ELSE IF (scheme==leapfrog_matsuno) THEN
87      CALL  leapfrog_matsuno_scheme
88    ELSE IF (scheme==adam_bashforth) THEN
89      CALL dissip(f_u,f_du)   
90      CALL adam_bashforth_scheme
91    ENDIF
92
93  ENDDO
94 
95  CONTAINS
96
97    SUBROUTINE Euler_scheme
98    IMPLICIT NONE
99      INTEGER :: ind
100     
101      DO ind=1,ndomain
102        h=f_h(ind) ; u=f_u(ind) ; 
103        dh=f_dh(ind) ; du=f_du(ind) ;
104        h(:)=h(:)+dt*dh(:)
105        u(:)=u(:)+dt*du(:)
106      ENDDO
107     
108    END SUBROUTINE Euler_scheme
109   
110    SUBROUTINE leapfrog_scheme
111    IMPLICIT NONE
112      INTEGER :: ind
113
114      DO ind=1,ndomain
115        h=f_h(ind)   ; u=f_u(ind)   ; 
116        hm1=f_hm1(ind) ; um1=f_um1(ind) ; 
117        hm2=f_hm2(ind) ; um2=f_um2(ind) ; 
118        dh=f_dh(ind) ; du=f_du(ind) ;
119
120        IF (it==0) THEN
121          hm2(:)=h(:) ; um2(:)=u(:) 
122
123          h(:)=h(:)+dt*dh(:)
124          u(:)=u(:)+dt*du(:)
125          hm1(:)=h(:) ;  um1(:)=u(:) 
126
127        ELSE
128       
129          h(:)=hm2(:)+2*dt*dh(:)
130          u(:)=um2(:)+2*dt*du(:)
131
132          hm2(:)=hm1(:) ; um2(:)=um1(:) 
133          hm1(:)=h(:) ; um1(:)=u(:) 
134        ENDIF
135         
136      ENDDO
137    END SUBROUTINE leapfrog_scheme 
138 
139    SUBROUTINE leapfrog_matsuno_scheme
140    IMPLICIT NONE
141      INTEGER :: ind
142
143      DO ind=1,ndomain
144        h=f_h(ind)   ; u=f_u(ind) 
145        hm1=f_hm1(ind) ; um1=f_um1(ind) 
146        hm2=f_hm2(ind) ; um2=f_um2(ind) 
147        dh=f_dh(ind) ; du=f_du(ind) ; 
148
149       
150        IF (MOD(it,matsuno_period+1)==0) THEN
151          hm1(:)=h(:) ; um1(:)=u(:) 
152          h(:)=hm1(:)+dt*dh(:)
153          u(:)=um1(:)+dt*du(:)
154
155        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
156
157          h(:)=hm1(:)+dt*dh(:)
158          u(:)=um1(:)+dt*du(:)
159
160          hm2(:)=hm1(:) ;  um2(:)=um1(:) 
161          hm1(:)=h(:) ;  um1(:)=u(:) 
162
163        ELSE
164
165          h(:)=hm2(:)+2*dt*dh(:)
166          u(:)=um2(:)+2*dt*du(:)
167
168          hm2(:)=hm1(:) ; um2(:)=um1(:) 
169          hm1(:)=h(:) ;  um1(:)=u(:) 
170
171        ENDIF
172     
173      ENDDO
174     
175    END SUBROUTINE leapfrog_matsuno_scheme 
176         
177    SUBROUTINE adam_bashforth_scheme
178    IMPLICIT NONE
179      INTEGER :: ind
180
181      DO ind=1,ndomain
182        h=f_h(ind)   ; u=f_u(ind)   
183        dh=f_dh(ind) ; du=f_du(ind) 
184        dhm1=f_dhm1(ind) ; dum1=f_dum1(ind) 
185        dhm2=f_dhm2(ind) ; dum2=f_dum2(ind) 
186
187        IF (it==0) THEN
188          dhm1(:)=dh(:) ; dum1(:)=du(:) 
189          dhm2(:)=dhm1(:) ; dum2(:)=dum1(:)
190        ENDIF
191             
192        h(:)=h(:)+dt*(23*dh(:)-16*dhm1(:)+5*dhm2(:))/12
193        u(:)=u(:)+dt*(23*du(:)-16*dum1(:)+5*dum2(:))/12
194
195        dhm2(:)=dhm1(:) ; dum2(:)=dum1(:) 
196        dhm1(:)=dh(:) ; dum1(:)=du(:) 
197
198      ENDDO
199     
200     
201    END SUBROUTINE adam_bashforth_scheme
202
203  END SUBROUTINE timeloop
204   
205   
206 
207END MODULE timeloop_sw_mod
Note: See TracBrowser for help on using the repository browser.