source: codes/icosagcm/trunk/src/timeloop_gcm.f90 @ 21

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

correction for compiling with gfortran (line too long)
improvement for splitting domain
Call twice transfert request for field u is no longer necessary

YM
File size: 8.0 KB
Line 
1MODULE timeloop_gcm_mod
2
3 
4CONTAINS
5 
6  SUBROUTINE timeloop
7  USE icosa
8  USE dissip_gcm_mod
9  USE caldyn_mod
10  USE theta2theta_rhodz_mod
11  USE etat0_mod
12  USE guided_mod
13  USE advect_tracer_mod
14 
15  IMPLICIT NONE
16  TYPE(t_field),POINTER :: f_phis(:)
17  TYPE(t_field),POINTER :: f_theta(:)
18  TYPE(t_field),POINTER :: f_q(:)
19  TYPE(t_field),POINTER :: f_dtheta(:)
20  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
21  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
22  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
23  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
24  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
25  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
26
27  REAL(rstd),POINTER :: phis(:)
28  REAL(rstd),POINTER :: q(:,:,:)
29  REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
30  REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
31  REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
32  REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
33  REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
34  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
35  REAL(rstd) :: dt
36  INTEGER :: ind
37  INTEGER :: it,i,j,n
38  CHARACTER(len=255) :: scheme
39  INTEGER :: matsuno_period
40  INTEGER :: itaumax
41  INTEGER :: write_period 
42  INTEGER :: itau_out
43
44  dt=90.
45  CALL getin('dt',dt)
46
47  itaumax=100
48  CALL getin('itaumax',itaumax)
49
50  write_period=0
51  CALL getin('write_period',write_period)
52  itau_out=INT(write_period/dt)
53 
54  scheme='adam_bashforth'
55  CALL getin('scheme',scheme)
56 
57  matsuno_period=5
58  CALL getin('matsuno_period',matsuno_period)
59  IF (TRIM(scheme)=='leapfrog') matsuno_period=itaumax+1
60
61  CALL allocate_field(f_phis,field_t,type_real)
62 
63  CALL allocate_field(f_ps,field_t,type_real)
64  CALL allocate_field(f_psm1,field_t,type_real)
65  CALL allocate_field(f_psm2,field_t,type_real)
66  CALL allocate_field(f_dps,field_t,type_real)
67  CALL allocate_field(f_dpsm1,field_t,type_real)
68  CALL allocate_field(f_dpsm2,field_t,type_real)
69
70  CALL allocate_field(f_u,field_u,type_real,llm)
71  CALL allocate_field(f_um1,field_u,type_real,llm)
72  CALL allocate_field(f_um2,field_u,type_real,llm)
73  CALL allocate_field(f_du,field_u,type_real,llm)
74  CALL allocate_field(f_dum1,field_u,type_real,llm)
75  CALL allocate_field(f_dum2,field_u,type_real,llm)
76
77  CALL allocate_field(f_theta,field_t,type_real,llm)
78  CALL allocate_field(f_dtheta,field_t,type_real,llm)
79
80  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
81
82  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
83  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
84  CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
85  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
86  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
87  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
88
89  CALL init_dissip(dt)
90  CALL init_caldyn(dt)
91  CALL init_guided(dt)
92!  CALL init_advect_tracer(dt)
93 
94  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
95 
96  DO it=0,itaumax
97    PRINT *,"It No :",It
98
99    CALL guided(it,f_ps,f_theta_rhodz,f_u,f_q)
100    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du)
101!    CALL advect_tracer(f_ps,f_u,f_q)
102   
103    SELECT CASE (TRIM(scheme))
104      CASE('euler')
105        CALL  euler_scheme
106
107      CASE ('leapfrog')
108        CALL leapfrog_scheme
109
110      CASE ('leapfrog_matsuno')
111        CALL  leapfrog_matsuno_scheme
112
113      CASE ('adam_bashforth')
114        CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz)
115        CALL adam_bashforth_scheme
116
117      CASE default
118        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
119               ' > options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>'
120        STOP
121       
122    END SELECT
123
124   
125    IF ( itau_out>0 .AND. MOD(it,itau_out)==0) THEN
126      CALL writefield("q",f_q)
127      CALL writefield("ps",f_ps)
128    ENDIF
129
130  ENDDO
131 
132  CONTAINS
133
134    SUBROUTINE Euler_scheme
135    IMPLICIT NONE
136      INTEGER :: ind
137     
138      DO ind=1,ndomain
139        ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
140        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
141        ps(:)=ps(:)+dt*dps(:)
142        u(:,:)=u(:,:)+dt*du(:,:)
143        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
144      ENDDO
145     
146    END SUBROUTINE Euler_scheme
147   
148    SUBROUTINE leapfrog_scheme
149    IMPLICIT NONE
150      INTEGER :: ind
151
152      DO ind=1,ndomain
153        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
154        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
155        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
156        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
157
158        IF (it==0) THEN
159          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
160
161          ps(:)=ps(:)+dt*dps(:)
162          u(:,:)=u(:,:)+dt*du(:,:)
163          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
164
165          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
166        ELSE
167       
168          ps(:)=psm2(:)+2*dt*dps(:)
169          u(:,:)=um2(:,:)+2*dt*du(:,:)
170          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
171
172          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
173          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
174        ENDIF
175         
176      ENDDO
177    END SUBROUTINE leapfrog_scheme 
178 
179    SUBROUTINE leapfrog_matsuno_scheme
180    IMPLICIT NONE
181      INTEGER :: ind
182
183      DO ind=1,ndomain
184        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
185        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
186        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
187        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
188
189       
190        IF (MOD(it,matsuno_period+1)==0) THEN
191          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
192          ps(:)=psm1(:)+dt*dps(:)
193          u(:,:)=um1(:,:)+dt*du(:,:)
194          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)
195
196        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
197
198          ps(:)=psm1(:)+dt*dps(:)
199          u(:,:)=um1(:,:)+dt*du(:,:)
200          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)
201
202          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
203          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
204
205        ELSE
206
207          ps(:)=psm2(:)+2*dt*dps(:)
208          u(:,:)=um2(:,:)+2*dt*du(:,:)
209          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
210
211          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
212          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
213
214        ENDIF
215     
216      ENDDO
217     
218    END SUBROUTINE leapfrog_matsuno_scheme 
219         
220    SUBROUTINE adam_bashforth_scheme
221    IMPLICIT NONE
222      INTEGER :: ind
223
224      DO ind=1,ndomain
225        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
226        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
227        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
228        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
229
230        IF (it==0) THEN
231          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
232          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
233        ENDIF
234             
235        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
236        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
237        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
238
239        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
240        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
241
242      ENDDO
243     
244     
245    END SUBROUTINE adam_bashforth_scheme
246
247  END SUBROUTINE timeloop
248   
249   
250 
251END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.