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

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

caldyn_adv is back
File size: 14.7 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  LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
49!  INTEGER :: itaumax
50!  REAL(rstd) ::write_period
51!  INTEGER    :: itau_out
52 
53!  dt=90.
54!  CALL getin('dt',dt)
55 
56!  itaumax=100
57!  CALL getin('itaumax',itaumax)
58
59!  run_length=dt*itaumax
60!  CALL getin('run_length',run_length)
61!  itaumax=run_length/dt
62!  PRINT *,'itaumax=',itaumax
63!  dt=dt/scale_factor
64
65! Trends
66  CALL allocate_field(f_dps,field_t,type_real)
67  CALL allocate_field(f_du,field_u,type_real,llm)
68  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
69! Model state at current time step (RK/MLF/Euler)
70  CALL allocate_field(f_phis,field_t,type_real)
71  CALL allocate_field(f_ps,field_t,type_real)
72  CALL allocate_field(f_u,field_u,type_real,llm)
73  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
74! Model state at previous time step (RK/MLF)
75  CALL allocate_field(f_psm1,field_t,type_real)
76  CALL allocate_field(f_um1,field_u,type_real,llm)
77  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
78! Tracers
79  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
80  CALL allocate_field(f_rhodz,field_t,type_real,llm)
81! Mass fluxes
82  CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
83  CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
84  CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
85  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
86
87  scheme_name='runge_kutta'
88  CALL getin('scheme',scheme_name)
89
90  SELECT CASE (TRIM(scheme_name))
91  CASE('euler')
92     scheme=euler
93     nb_stage=1
94  CASE ('runge_kutta')
95     scheme=rk4
96     nb_stage=4
97  CASE ('leapfrog_matsuno')
98     scheme=mlf
99     matsuno_period=5
100     CALL getin('matsuno_period',matsuno_period)
101     nb_stage=matsuno_period+1
102     ! Model state 2 time steps ago (MLF)
103     CALL allocate_field(f_psm2,field_t,type_real)
104     CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
105     CALL allocate_field(f_um2,field_u,type_real,llm)
106  CASE default
107     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
108          ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
109     STOP
110  END SELECT
111 
112!  write_period=0
113!  CALL getin('write_period',write_period)
114!  write_period=write_period/scale_factor
115!  itau_out=FLOOR(.5+write_period/dt)
116!  PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
117 
118! Trends at previous time steps needed only by Adams-Bashforth
119!  CALL allocate_field(f_dpsm1,field_t,type_real)
120!  CALL allocate_field(f_dpsm2,field_t,type_real)
121!  CALL allocate_field(f_dum1,field_u,type_real,llm)
122!  CALL allocate_field(f_dum2,field_u,type_real,llm)
123!  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
124!  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
125
126!  CALL allocate_field(f_theta,field_t,type_real,llm)
127!  CALL allocate_field(f_dtheta,field_t,type_real,llm)
128
129  CALL init_dissip
130  CALL init_caldyn
131  CALL init_guided
132  CALL init_advect_tracer
133  CALL init_physics
134 
135  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
136  CALL writefield("phis",f_phis,once=.TRUE.)
137  CALL transfert_request(f_q,req_i1) 
138
139  DO ind=1,ndomain
140     CALL swap_dimensions(ind)
141     CALL swap_geometry(ind)
142     rhodz=f_rhodz(ind); ps=f_ps(ind)
143     CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
144  END DO
145  fluxt_zero=.TRUE.
146
147  ! check that rhodz is consistent with ps
148  CALL transfert_request(f_rhodz,req_i1)
149  CALL transfert_request(f_ps,req_i1)
150  DO ind=1,ndomain
151     CALL swap_dimensions(ind)
152     CALL swap_geometry(ind)
153     rhodz=f_rhodz(ind); ps=f_ps(ind)
154     CALL compute_rhodz(.FALSE., ps, rhodz)   
155  END DO
156 
157  DO it=0,itaumax
158
159    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
160    IF (mod(it,itau_out)==0 ) THEN
161      CALL writefield("q",f_q)
162      CALL update_time_counter(dt*it)
163    ENDIF
164   
165    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
166
167    DO stage=1,nb_stage
168       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
169            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
170            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
171       SELECT CASE (scheme)
172       CASE(euler)
173          CALL euler_scheme(.TRUE.)
174       CASE (rk4)
175          CALL rk_scheme(stage)
176       CASE (mlf)
177          CALL  leapfrog_matsuno_scheme(stage)
178         
179          !      CASE ('leapfrog')
180          !        CALL leapfrog_scheme
181          !
182          !      CASE ('adam_bashforth')
183          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
184          !        CALL adam_bashforth_scheme
185       CASE DEFAULT
186          STOP
187       END SELECT
188    END DO
189
190    CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
191    CALL euler_scheme(.FALSE.)
192
193    IF(MOD(it+1,itau_adv)==0) THEN
194       CALL transfert_request(f_wfluxt,req_i1) ! FIXME
195!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
196
197       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
198       fluxt_zero=.TRUE.
199
200       ! FIXME : check that rhodz is consistent with ps
201       CALL transfert_request(f_rhodz,req_i1)
202       CALL transfert_request(f_ps,req_i1)
203       CALL transfert_request(f_dps,req_i1) ! FIXME
204       CALL transfert_request(f_wflux,req_i1) ! FIXME
205       DO ind=1,ndomain
206          CALL swap_dimensions(ind)
207          CALL swap_geometry(ind)
208          rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 
209          wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
210          CALL compute_rhodz(.FALSE., ps, rhodz)   
211       END DO
212
213    END IF
214
215!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
216    ENDDO
217 
218  CONTAINS
219
220    SUBROUTINE Euler_scheme(with_dps)
221    IMPLICIT NONE
222    LOGICAL :: with_dps
223    INTEGER :: ind
224    DO ind=1,ndomain
225       CALL swap_dimensions(ind)
226       CALL swap_geometry(ind)
227       IF(with_dps) THEN
228          ps=f_ps(ind) ; dps=f_dps(ind) ; 
229          ps(:)=ps(:)+dt*dps(:)
230          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
231          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
232          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
233       END IF
234       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
235       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
236       u(:,:)=u(:,:)+dt*du(:,:)
237       theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
238    ENDDO
239
240    END SUBROUTINE Euler_scheme
241
242    SUBROUTINE RK_scheme(stage)
243      IMPLICIT NONE
244      INTEGER :: ind, stage
245      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
246      REAL(rstd) :: tau
247
248      tau = dt*coef(stage)
249     
250      DO ind=1,ndomain
251         CALL swap_dimensions(ind)
252         CALL swap_geometry(ind)
253         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
254         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
255         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
256         
257         IF (stage==1) THEN ! first stage : save model state in XXm1
258            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
259         END IF
260         ! updates are of the form x1 := x0 + tau*f(x1)
261         ps(:)=psm1(:)+tau*dps(:)
262         u(:,:)=um1(:,:)+tau*du(:,:)
263         theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
264         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
265            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
266            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
267            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
268         END IF
269      END DO
270    END SUBROUTINE RK_scheme
271
272    SUBROUTINE leapfrog_scheme
273    IMPLICIT NONE
274      INTEGER :: ind
275
276      DO ind=1,ndomain
277        CALL swap_dimensions(ind)
278        CALL swap_geometry(ind)
279        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
280        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
281        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
282        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
283
284        IF (it==0) THEN
285          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
286
287          ps(:)=ps(:)+dt*dps(:)
288          u(:,:)=u(:,:)+dt*du(:,:)
289          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
290
291          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
292        ELSE
293       
294          ps(:)=psm2(:)+2*dt*dps(:)
295          u(:,:)=um2(:,:)+2*dt*du(:,:)
296          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
297
298          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
299          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
300        ENDIF
301         
302      ENDDO
303    END SUBROUTINE leapfrog_scheme 
304 
305    SUBROUTINE leapfrog_matsuno_scheme(stage)
306    IMPLICIT NONE
307    INTEGER :: ind, stage
308    REAL :: tau
309    tau = dt/nb_stage
310      DO ind=1,ndomain
311        CALL swap_dimensions(ind)
312        CALL swap_geometry(ind)
313
314        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
315        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
316        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
317        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
318
319       
320!        IF (MOD(it,matsuno_period+1)==0) THEN
321        IF (stage==1) THEN
322          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
323          ps(:)=psm1(:)+tau*dps(:)
324          u(:,:)=um1(:,:)+tau*du(:,:)
325          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
326
327!        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
328        ELSE IF (stage==2) THEN
329
330          ps(:)=psm1(:)+tau*dps(:)
331          u(:,:)=um1(:,:)+tau*du(:,:)
332          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
333
334          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
335          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
336
337        ELSE
338
339          ps(:)=psm2(:)+2*tau*dps(:)
340          u(:,:)=um2(:,:)+2*tau*du(:,:)
341          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
342
343          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
344          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
345
346        ENDIF
347     
348      ENDDO
349     
350    END SUBROUTINE leapfrog_matsuno_scheme 
351         
352    SUBROUTINE adam_bashforth_scheme
353    IMPLICIT NONE
354      INTEGER :: ind
355
356      DO ind=1,ndomain
357        CALL swap_dimensions(ind)
358        CALL swap_geometry(ind)
359        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
360        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
361        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
362        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
363
364        IF (it==0) THEN
365          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
366          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
367        ENDIF
368             
369        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
370        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
371        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
372
373        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
374        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
375
376      ENDDO     
377     
378    END SUBROUTINE adam_bashforth_scheme
379
380  END SUBROUTINE timeloop   
381
382  SUBROUTINE compute_rhodz(comp, ps, rhodz)
383    USE icosa
384    USE disvert_mod
385    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
386    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
387    REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
388    REAL(rstd) :: m, err
389    INTEGER :: l,i,j,ij,dd
390    err=0.
391    IF(comp) THEN
392       dd=1
393    ELSE
394!       dd=-1
395       dd=0
396    END IF
397
398    DO l = 1, llm
399       DO j=jj_begin-dd,jj_end+dd
400          DO i=ii_begin-dd,ii_end+dd
401             ij=(j-1)*iim+i
402             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
403             IF(comp) THEN
404                rhodz(ij,l) = m
405             ELSE
406                err = MAX(err,abs(m-rhodz(ij,l)))
407             END IF
408          ENDDO
409       ENDDO
410    ENDDO
411
412    IF(.NOT. comp) THEN
413       IF(err>1e-10) THEN
414          PRINT *, 'Discrepancy between ps and rhodz detected', err
415          STOP
416       ELSE
417!          PRINT *, 'No discrepancy between ps and rhodz detected'
418       END IF
419    END IF
420
421  END SUBROUTINE compute_rhodz
422
423  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
424    USE icosa
425    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
426    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
427    REAL(rstd), INTENT(IN) :: tau
428    LOGICAL, INTENT(INOUT) :: fluxt_zero
429    IF(fluxt_zero) THEN
430!       PRINT *, 'Accumulating fluxes (first)'
431       fluxt_zero=.FALSE.
432       hfluxt = tau*hflux
433       wfluxt = tau*wflux
434    ELSE
435!       PRINT *, 'Accumulating fluxes (next)'
436       hfluxt = hfluxt + tau*hflux
437       wfluxt = wfluxt + tau*wflux
438    END IF
439  END SUBROUTINE accumulate_fluxes
440 
441END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.