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

Last change on this file since 133 was 133, checked in by dubos, 11 years ago

Working on dyn/adv coupling
File size: 12.8 KB
Line 
1MODULE timeloop_gcm_mod
2
3  PRIVATE
4 
5  PUBLIC :: timeloop
6
7  INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
8
9CONTAINS
10 
11  SUBROUTINE timeloop
12  USE icosa
13  USE dissip_gcm_mod
14  USE caldyn_mod
15  USE etat0_mod
16  USE guided_mod
17  USE advect_tracer_mod
18  USE physics_mod
19  USE mpipara
20 
21  IMPLICIT NONE
22  TYPE(t_field),POINTER :: f_phis(:)
23!  TYPE(t_field),POINTER :: f_theta(:)
24  TYPE(t_field),POINTER :: f_q(:)
25  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
26  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
27  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
28  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
29  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
30  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
31  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
32  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
33
34  REAL(rstd),POINTER :: phis(:)
35  REAL(rstd),POINTER :: q(:,:,:)
36  REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
37  REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
38  REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
39  REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
40  REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
41  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
42  REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
43
44!  REAL(rstd) :: dt, run_length
45  INTEGER :: ind
46  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
47  CHARACTER(len=255) :: scheme_name
48!  INTEGER :: itaumax
49!  REAL(rstd) ::write_period
50!  INTEGER    :: itau_out
51 
52!  dt=90.
53!  CALL getin('dt',dt)
54 
55!  itaumax=100
56!  CALL getin('itaumax',itaumax)
57
58!  run_length=dt*itaumax
59!  CALL getin('run_length',run_length)
60!  itaumax=run_length/dt
61!  PRINT *,'itaumax=',itaumax
62!  dt=dt/scale_factor
63
64! Trends
65  CALL allocate_field(f_dps,field_t,type_real)
66  CALL allocate_field(f_du,field_u,type_real,llm)
67  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
68! Model state at current time step (RK/MLF/Euler)
69  CALL allocate_field(f_phis,field_t,type_real)
70  CALL allocate_field(f_ps,field_t,type_real)
71  CALL allocate_field(f_u,field_u,type_real,llm)
72  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
73! Model state at previous time step (RK/MLF)
74  CALL allocate_field(f_psm1,field_t,type_real)
75  CALL allocate_field(f_um1,field_u,type_real,llm)
76  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
77! Tracers
78  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
79  CALL allocate_field(f_rhodz,field_t,type_real,llm)
80  CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
81  CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
82  CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
83  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
84
85  scheme_name='runge_kutta'
86  CALL getin('scheme',scheme_name)
87
88  SELECT CASE (TRIM(scheme_name))
89  CASE('euler')
90     scheme=euler
91     nb_stage=1
92  CASE ('runge_kutta')
93     scheme=rk4
94     nb_stage=4
95  CASE ('leapfrog_matsuno')
96     scheme=mlf
97     matsuno_period=5
98     CALL getin('matsuno_period',matsuno_period)
99     nb_stage=matsuno_period+1
100     ! Model state 2 time steps ago (MLF)
101     CALL allocate_field(f_psm2,field_t,type_real)
102     CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
103     CALL allocate_field(f_um2,field_u,type_real,llm)
104  CASE default
105     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
106          ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
107     STOP
108  END SELECT
109 
110!  write_period=0
111!  CALL getin('write_period',write_period)
112!  write_period=write_period/scale_factor
113!  itau_out=FLOOR(.5+write_period/dt)
114!  PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
115 
116! Trends at previous time steps needed only by Adams-Bashforth
117!  CALL allocate_field(f_dpsm1,field_t,type_real)
118!  CALL allocate_field(f_dpsm2,field_t,type_real)
119!  CALL allocate_field(f_dum1,field_u,type_real,llm)
120!  CALL allocate_field(f_dum2,field_u,type_real,llm)
121!  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
122!  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
123
124!  CALL allocate_field(f_theta,field_t,type_real,llm)
125!  CALL allocate_field(f_dtheta,field_t,type_real,llm)
126
127  CALL init_dissip
128  CALL init_caldyn
129  CALL init_guided
130  CALL init_advect_tracer
131  CALL init_physics
132 
133  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
134  CALL writefield("phis",f_phis,once=.TRUE.)
135  CALL transfert_request(f_q,req_i1) 
136
137  DO ind=1,ndomain
138     CALL swap_dimensions(ind)
139     CALL swap_geometry(ind)
140     hfluxt=f_hfluxt(ind); hfluxt=0. ! clear accumulated fluxes
141     wfluxt=f_hfluxt(ind); wfluxt=0.
142     rhodz=f_rhodz(ind); ps=f_ps(ind)
143     CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
144  END DO
145
146  DO it=0,itaumax
147
148    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
149    IF (mod(it,itau_out)==0 ) THEN
150      CALL writefield("q",f_q)
151      CALL update_time_counter(dt*it)
152    ENDIF
153   
154!    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
155
156    DO stage=1,nb_stage
157!       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
158!            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
159!            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
160       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
161            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
162            f_dps, f_dtheta_rhodz, f_du)
163       SELECT CASE (scheme)
164       CASE(euler)
165          CALL euler_scheme(.TRUE.)
166       CASE (rk4)
167          CALL rk_scheme(stage)
168       CASE (mlf)
169          CALL  leapfrog_matsuno_scheme(stage)
170         
171          !      CASE ('leapfrog')
172          !        CALL leapfrog_scheme
173          !
174          !      CASE ('adam_bashforth')
175          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
176          !        CALL adam_bashforth_scheme
177       CASE DEFAULT
178          STOP
179       END SELECT
180    END DO
181
182    CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
183    CALL euler_scheme(.FALSE.)
184
185    IF(MOD(it+1,itau_adv)==0) THEN
186!       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
187       DO ind=1,ndomain
188          hfluxt=f_hfluxt(ind)
189          wfluxt=f_hfluxt(ind)
190          hfluxt=0.               ! clear accumulated fluxes
191          wfluxt=0.
192       END DO
193    END IF
194
195!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
196    ENDDO
197 
198  CONTAINS
199
200    SUBROUTINE Euler_scheme(with_dps)
201    IMPLICIT NONE
202    LOGICAL :: with_dps
203    INTEGER :: ind
204     
205    DO ind=1,ndomain
206       IF(with_dps) THEN
207          ps=f_ps(ind) ; dps=f_dps(ind) ; 
208          ps(:)=ps(:)+dt*dps(:)
209          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
210          wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
211          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt)
212       END IF
213       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
214       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
215       u(:,:)=u(:,:)+dt*du(:,:)
216       theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
217    ENDDO
218
219    END SUBROUTINE Euler_scheme
220
221    SUBROUTINE RK_scheme(stage)
222      IMPLICIT NONE
223      INTEGER :: ind, stage
224      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
225      REAL(rstd) :: tau
226
227      tau = dt*coef(stage)
228     
229      DO ind=1,ndomain
230         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
231         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
232         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
233         
234         IF (stage==1) THEN ! first stage : save model state in XXm1
235            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
236         END IF
237         ! updates are of the form x1 := x0 + tau*f(x1)
238         ps(:)=psm1(:)+tau*dps(:)
239         u(:,:)=um1(:,:)+tau*du(:,:)
240         theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
241         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
242            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
243            wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
244            CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt)
245         END IF
246      END DO
247    END SUBROUTINE RK_scheme
248
249    SUBROUTINE leapfrog_scheme
250    IMPLICIT NONE
251      INTEGER :: ind
252
253      DO ind=1,ndomain
254        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
255        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
256        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
257        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
258
259        IF (it==0) THEN
260          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
261
262          ps(:)=ps(:)+dt*dps(:)
263          u(:,:)=u(:,:)+dt*du(:,:)
264          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
265
266          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
267        ELSE
268       
269          ps(:)=psm2(:)+2*dt*dps(:)
270          u(:,:)=um2(:,:)+2*dt*du(:,:)
271          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
272
273          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
274          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
275        ENDIF
276         
277      ENDDO
278    END SUBROUTINE leapfrog_scheme 
279 
280    SUBROUTINE leapfrog_matsuno_scheme(stage)
281    IMPLICIT NONE
282    INTEGER :: ind, stage
283    REAL :: tau
284    tau = dt/nb_stage
285      DO ind=1,ndomain
286        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
287        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
288        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
289        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
290
291       
292!        IF (MOD(it,matsuno_period+1)==0) THEN
293        IF (stage==1) THEN
294          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
295          ps(:)=psm1(:)+tau*dps(:)
296          u(:,:)=um1(:,:)+tau*du(:,:)
297          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
298
299!        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
300        ELSE IF (stage==2) THEN
301
302          ps(:)=psm1(:)+tau*dps(:)
303          u(:,:)=um1(:,:)+tau*du(:,:)
304          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
305
306          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
307          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
308
309        ELSE
310
311          ps(:)=psm2(:)+2*tau*dps(:)
312          u(:,:)=um2(:,:)+2*tau*du(:,:)
313          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
314
315          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
316          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
317
318        ENDIF
319     
320      ENDDO
321     
322    END SUBROUTINE leapfrog_matsuno_scheme 
323         
324    SUBROUTINE adam_bashforth_scheme
325    IMPLICIT NONE
326      INTEGER :: ind
327
328      DO ind=1,ndomain
329        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
330        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
331        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
332        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
333
334        IF (it==0) THEN
335          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
336          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
337        ENDIF
338             
339        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
340        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
341        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
342
343        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
344        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
345
346      ENDDO     
347     
348    END SUBROUTINE adam_bashforth_scheme
349
350  END SUBROUTINE timeloop   
351
352  SUBROUTINE compute_rhodz(ps, rhodz)
353    USE icosa
354    USE disvert_mod
355    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
356    REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
357    INTEGER :: l,i,j,ij
358    DO l = 1, llm
359       DO j=jj_begin-1,jj_end+1
360          DO i=ii_begin-1,ii_end+1
361             ij=(j-1)*iim+i
362             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
363          ENDDO
364       ENDDO
365    ENDDO
366  END SUBROUTINE compute_rhodz
367
368  SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau)
369    USE icosa
370    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
371    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
372    REAL(rstd), INTENT(IN) :: tau
373    hfluxt = hfluxt + tau*hflux
374    wfluxt = wfluxt + tau*wflux
375  END SUBROUTINE accumulate_fluxes
376 
377END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.