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

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

Removed duplicate write_field
Solved compilation issues with cpp on Bluefire (use GNU cpp)
Tested : Test case 3 with nbp=20 and interpolation to lat-lon.
File size: 8.1 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, run_length
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  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  write_period=0
57  CALL getin('write_period',write_period)
58  write_period=write_period/scale_factor
59 
60  itau_out=INT(write_period/dt)
61 
62  scheme='adam_bashforth'
63  CALL getin('scheme',scheme)
64 
65  matsuno_period=5
66  CALL getin('matsuno_period',matsuno_period)
67  IF (TRIM(scheme)=='leapfrog') matsuno_period=itaumax+1
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(dt)
98  CALL init_caldyn(dt)
99  CALL init_guided(dt)
100  CALL init_advect_tracer(dt)
101 
102  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
103 
104  DO it=0,itaumax
105    PRINT *,"It No :",It,"   t :",dt*It
106
107    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
108    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du)
109    CALL advect_tracer(f_ps,f_u,f_q)
110   
111    SELECT CASE (TRIM(scheme))
112      CASE('euler')
113        CALL  euler_scheme
114
115      CASE ('leapfrog')
116        CALL leapfrog_scheme
117
118      CASE ('leapfrog_matsuno')
119        CALL  leapfrog_matsuno_scheme
120
121      CASE ('adam_bashforth')
122        CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz)
123        CALL adam_bashforth_scheme
124
125      CASE default
126        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
127               ' > options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>'
128        STOP
129       
130    END SELECT
131
132 ENDDO
133 
134  CONTAINS
135
136    SUBROUTINE Euler_scheme
137    IMPLICIT NONE
138      INTEGER :: ind
139     
140      DO ind=1,ndomain
141        ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
142        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
143        ps(:)=ps(:)+dt*dps(:)
144        u(:,:)=u(:,:)+dt*du(:,:)
145        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
146      ENDDO
147     
148    END SUBROUTINE Euler_scheme
149   
150    SUBROUTINE leapfrog_scheme
151    IMPLICIT NONE
152      INTEGER :: ind
153
154      DO ind=1,ndomain
155        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
156        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
157        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
158        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
159
160        IF (it==0) THEN
161          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
162
163          ps(:)=ps(:)+dt*dps(:)
164          u(:,:)=u(:,:)+dt*du(:,:)
165          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
166
167          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
168        ELSE
169       
170          ps(:)=psm2(:)+2*dt*dps(:)
171          u(:,:)=um2(:,:)+2*dt*du(:,:)
172          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
173
174          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
175          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
176        ENDIF
177         
178      ENDDO
179    END SUBROUTINE leapfrog_scheme 
180 
181    SUBROUTINE leapfrog_matsuno_scheme
182    IMPLICIT NONE
183      INTEGER :: ind
184
185      DO ind=1,ndomain
186        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
187        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
188        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
189        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
190
191       
192        IF (MOD(it,matsuno_period+1)==0) THEN
193          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
194          ps(:)=psm1(:)+dt*dps(:)
195          u(:,:)=um1(:,:)+dt*du(:,:)
196          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)
197
198        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
199
200          ps(:)=psm1(:)+dt*dps(:)
201          u(:,:)=um1(:,:)+dt*du(:,:)
202          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)
203
204          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
205          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
206
207        ELSE
208
209          ps(:)=psm2(:)+2*dt*dps(:)
210          u(:,:)=um2(:,:)+2*dt*du(:,:)
211          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
212
213          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
214          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
215
216        ENDIF
217     
218      ENDDO
219     
220    END SUBROUTINE leapfrog_matsuno_scheme 
221         
222    SUBROUTINE adam_bashforth_scheme
223    IMPLICIT NONE
224      INTEGER :: ind
225
226      DO ind=1,ndomain
227        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
228        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
229        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
230        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
231
232        IF (it==0) THEN
233          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
234          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
235        ENDIF
236             
237        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
238        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
239        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
240
241        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
242        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
243
244      ENDDO
245     
246     
247    END SUBROUTINE adam_bashforth_scheme
248
249  END SUBROUTINE timeloop
250   
251   
252 
253END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.