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

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

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