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

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

Put time variable : dt, itaumax, write_period, itau_out in the time module

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