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

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

Correction for dcmip moist physics
Temperature is ouput instead of virtual temperature.

YM
File size: 8.5 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 ('leapfrog_matsuno')
127        CALL  leapfrog_matsuno_scheme
128
129      CASE ('adam_bashforth')
130        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
131        CALL adam_bashforth_scheme
132
133      CASE default
134        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
135               ' > options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>'
136        STOP
137       
138    END SELECT
139
140    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
141   
142 ENDDO
143 
144  CONTAINS
145
146    SUBROUTINE Euler_scheme
147    IMPLICIT NONE
148      INTEGER :: ind
149     
150      DO ind=1,ndomain
151        ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
152        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
153        ps(:)=ps(:)+dt*dps(:)
154        u(:,:)=u(:,:)+dt*du(:,:)
155        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
156      ENDDO
157     
158    END SUBROUTINE Euler_scheme
159   
160    SUBROUTINE leapfrog_scheme
161    IMPLICIT NONE
162      INTEGER :: ind
163
164      DO ind=1,ndomain
165        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
166        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
167        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
168        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
169
170        IF (it==0) THEN
171          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
172
173          ps(:)=ps(:)+dt*dps(:)
174          u(:,:)=u(:,:)+dt*du(:,:)
175          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
176
177          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
178        ELSE
179       
180          ps(:)=psm2(:)+2*dt*dps(:)
181          u(:,:)=um2(:,:)+2*dt*du(:,:)
182          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
183
184          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
185          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
186        ENDIF
187         
188      ENDDO
189    END SUBROUTINE leapfrog_scheme 
190 
191    SUBROUTINE leapfrog_matsuno_scheme
192    IMPLICIT NONE
193      INTEGER :: ind
194
195      DO ind=1,ndomain
196        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
197        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
198        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
199        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
200
201       
202        IF (MOD(it,matsuno_period+1)==0) THEN
203          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
204          ps(:)=psm1(:)+dt*dps(:)
205          u(:,:)=um1(:,:)+dt*du(:,:)
206          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)
207
208        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
209
210          ps(:)=psm1(:)+dt*dps(:)
211          u(:,:)=um1(:,:)+dt*du(:,:)
212          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)
213
214          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
215          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
216
217        ELSE
218
219          ps(:)=psm2(:)+2*dt*dps(:)
220          u(:,:)=um2(:,:)+2*dt*du(:,:)
221          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
222
223          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
224          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
225
226        ENDIF
227     
228      ENDDO
229     
230    END SUBROUTINE leapfrog_matsuno_scheme 
231         
232    SUBROUTINE adam_bashforth_scheme
233    IMPLICIT NONE
234      INTEGER :: ind
235
236      DO ind=1,ndomain
237        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
238        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
239        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
240        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
241
242        IF (it==0) THEN
243          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
244          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
245        ENDIF
246             
247        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
248        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
249        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
250
251        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
252        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
253
254      ENDDO     
255     
256    END SUBROUTINE adam_bashforth_scheme
257
258  END SUBROUTINE timeloop   
259 
260END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.