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
RevLine 
[12]1MODULE timeloop_gcm_mod
2 
3CONTAINS
4 
5  SUBROUTINE timeloop
[19]6  USE icosa
[15]7  USE dissip_gcm_mod
[17]8  USE caldyn_mod
[12]9  USE etat0_mod
[17]10  USE guided_mod
11  USE advect_tracer_mod
[81]12  USE physics_mod
[17]13 
[12]14  IMPLICIT NONE
15  TYPE(t_field),POINTER :: f_phis(:)
16  TYPE(t_field),POINTER :: f_theta(:)
[17]17  TYPE(t_field),POINTER :: f_q(:)
[15]18  TYPE(t_field),POINTER :: f_dtheta(:)
[12]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(:)
[17]27  REAL(rstd),POINTER :: q(:,:,:)
[12]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(:,:)
[50]34
[98]35!  REAL(rstd) :: dt, run_length
[12]36  INTEGER :: ind
37  INTEGER :: it,i,j,n
[17]38  CHARACTER(len=255) :: scheme
[12]39  INTEGER :: matsuno_period
[98]40!  INTEGER :: itaumax
41!  REAL(rstd) ::write_period
42!  INTEGER    :: itau_out
[63]43 
[98]44!  dt=90.
45!  CALL getin('dt',dt)
[32]46 
[98]47!  itaumax=100
48!  CALL getin('itaumax',itaumax)
[17]49
[98]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
[45]55
[17]56  scheme='adam_bashforth'
[12]57  CALL getin('scheme',scheme)
58 
59  matsuno_period=5
60  CALL getin('matsuno_period',matsuno_period)
[17]61  IF (TRIM(scheme)=='leapfrog') matsuno_period=itaumax+1
[12]62
[98]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
[63]68 
[12]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
[15]85  CALL allocate_field(f_theta,field_t,type_real,llm)
86  CALL allocate_field(f_dtheta,field_t,type_real,llm)
87
[17]88  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
89
[12]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
[98]97  CALL init_dissip
98  CALL init_caldyn
99  CALL init_guided
100  CALL init_advect_tracer
101  CALL init_physics
[12]102 
[17]103  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
[97]104  CALL writefield("phis",f_phis,once=.TRUE.)
[81]105  CALL transfert_request(f_q,req_i1) 
[50]106
[12]107  DO it=0,itaumax
[81]108
[45]109    PRINT *,"It No :",It,"   t :",dt*It
[63]110    IF (mod(it,itau_out)==0 ) THEN
111      CALL writefield("q",f_q)
[81]112      CALL update_time_counter(dt*it)
[63]113    ENDIF
114   
[25]115    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
[110]116    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
[25]117    CALL advect_tracer(f_ps,f_u,f_q)
[17]118   
119    SELECT CASE (TRIM(scheme))
120      CASE('euler')
121        CALL  euler_scheme
[12]122
[17]123      CASE ('leapfrog')
124        CALL leapfrog_scheme
125
126      CASE ('leapfrog_matsuno')
127        CALL  leapfrog_matsuno_scheme
128
129      CASE ('adam_bashforth')
[109]130        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
[17]131        CALL adam_bashforth_scheme
132
133      CASE default
[21]134        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
135               ' > options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>'
[17]136        STOP
137       
138    END SELECT
[110]139
140    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
[81]141   
[46]142 ENDDO
[12]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
[50]254      ENDDO     
[12]255     
256    END SUBROUTINE adam_bashforth_scheme
257
[50]258  END SUBROUTINE timeloop   
[12]259 
260END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.