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

Last change on this file since 124 was 124, checked in by dubos, 11 years ago

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