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

Last change on this file since 25 was 25, checked in by dubos, 12 years ago

Minor changes :
caldyn_sw.f90, advect_tracer.f90
icosa_mod.f90 : added parameters for NCAR test cases needing global scope
guided_mod.f90 : CALL to guided_ncar now takes tt=it*dt instead of it as input

Significant changes :
timeloop_gcm.f90 : re-activated CALL to advection scheme
disvert_ncar.f90,
etat0_ncar.f90
guided_ncar_mod.f90 : simplification, introduced several getin(...), update due to recent changes in advection test cases (deformational flow, Hadley cell)
run_adv.def : new keys, reorganized for legibility

Tests :
icosa_gcm.exe tested with ncar_adv_shape=const and ncar_adv_wind=solid,deform,hadley.
q1=1 maintained to machine accuracy. Surface pressure slightly oscillates as expected.

FIXME : Tests by Sarvesh with revision 24 show incorrect advection of cosine bell by solid-body rotation. Not fixed.
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*dt,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.