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

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

Merge advection scheme from sarvesh in standard version

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