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

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

Simplify the management of the module.

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