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

Last change on this file since 146 was 146, checked in by ymipsl, 11 years ago

Set constant sign for wind way :
ne(ij,right)==ne_right=1
ne(ij,rup)==ne_rup=-1
ne(ij,lup)==ne_lup=1
ne(ij,left)==ne_left=-1
ne(ij,ldown)==ne_ldown=1
ne(ij,rdown)==ne_rdown=-1

Modified transfert function to be compliant for this convention.

YM
File size: 15.1 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  USE trace
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
225    CALL trace_start("Euler_scheme") 
226
227    DO ind=1,ndomain
228       CALL swap_dimensions(ind)
229       CALL swap_geometry(ind)
230       IF(with_dps) THEN
231          ps=f_ps(ind) ; dps=f_dps(ind) ; 
232          ps(:)=ps(:)+dt*dps(:)
233          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
234          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
235          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
236       END IF
237       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
238       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
239       u(:,:)=u(:,:)+dt*du(:,:)
240       theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
241    ENDDO
242
243    CALL trace_end("Euler_scheme") 
244
245    END SUBROUTINE Euler_scheme
246
247    SUBROUTINE RK_scheme(stage)
248      IMPLICIT NONE
249      INTEGER :: ind, stage
250      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
251      REAL(rstd) :: tau
252 
253      CALL trace_start("RK_scheme") 
254
255      tau = dt*coef(stage)
256     
257      DO ind=1,ndomain
258         CALL swap_dimensions(ind)
259         CALL swap_geometry(ind)
260         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
261         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
262         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
263         
264         IF (stage==1) THEN ! first stage : save model state in XXm1
265            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
266         END IF
267         ! updates are of the form x1 := x0 + tau*f(x1)
268         ps(:)=psm1(:)+tau*dps(:)
269         u(:,:)=um1(:,:)+tau*du(:,:)
270         theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
271         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
272            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
273            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
274            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
275         END IF
276      END DO
277     
278      CALL trace_end("RK_scheme")
279     
280    END SUBROUTINE RK_scheme
281
282    SUBROUTINE leapfrog_scheme
283    IMPLICIT NONE
284      INTEGER :: ind
285
286      CALL trace_start("leapfrog_scheme")
287       
288      DO ind=1,ndomain
289        CALL swap_dimensions(ind)
290        CALL swap_geometry(ind)
291        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
292        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
293        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
294        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
295
296        IF (it==0) THEN
297          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 
298
299          ps(:)=ps(:)+dt*dps(:)
300          u(:,:)=u(:,:)+dt*du(:,:)
301          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
302
303          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
304        ELSE
305       
306          ps(:)=psm2(:)+2*dt*dps(:)
307          u(:,:)=um2(:,:)+2*dt*du(:,:)
308          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)
309
310          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
311          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
312        ENDIF
313         
314      ENDDO
315
316      CALL trace_end("leapfrog_scheme")
317
318    END SUBROUTINE leapfrog_scheme 
319 
320    SUBROUTINE leapfrog_matsuno_scheme(stage)
321    IMPLICIT NONE
322    INTEGER :: ind, stage
323    REAL :: tau
324
325      CALL trace_start("leapfrog_matsuno_scheme")
326   
327      tau = dt/nb_stage
328      DO ind=1,ndomain
329        CALL swap_dimensions(ind)
330        CALL swap_geometry(ind)
331
332        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
333        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
334        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
335        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
336
337       
338!        IF (MOD(it,matsuno_period+1)==0) THEN
339        IF (stage==1) THEN
340          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
341          ps(:)=psm1(:)+tau*dps(:)
342          u(:,:)=um1(:,:)+tau*du(:,:)
343          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
344
345!        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
346        ELSE IF (stage==2) THEN
347
348          ps(:)=psm1(:)+tau*dps(:)
349          u(:,:)=um1(:,:)+tau*du(:,:)
350          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
351
352          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
353          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
354
355        ELSE
356
357          ps(:)=psm2(:)+2*tau*dps(:)
358          u(:,:)=um2(:,:)+2*tau*du(:,:)
359          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
360
361          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
362          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
363
364        ENDIF
365     
366      ENDDO
367      CALL trace_end("leapfrog_matsuno_scheme")
368     
369    END SUBROUTINE leapfrog_matsuno_scheme 
370         
371    SUBROUTINE adam_bashforth_scheme
372    IMPLICIT NONE
373      INTEGER :: ind
374
375      CALL trace_start("adam_bashforth_scheme")
376
377      DO ind=1,ndomain
378        CALL swap_dimensions(ind)
379        CALL swap_geometry(ind)
380        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
381        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
382        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
383        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)
384
385        IF (it==0) THEN
386          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
387          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
388        ENDIF
389             
390        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
391        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
392        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12
393
394        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
395        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
396
397      ENDDO     
398
399      CALL trace_end("adam_bashforth_scheme")
400     
401    END SUBROUTINE adam_bashforth_scheme
402
403  END SUBROUTINE timeloop   
404
405  SUBROUTINE compute_rhodz(comp, ps, rhodz)
406    USE icosa
407    USE disvert_mod
408    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
409    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
410    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
411    REAL(rstd) :: m, err
412    INTEGER :: l,i,j,ij,dd
413    err=0.
414    IF(comp) THEN
415       dd=1
416    ELSE
417!       dd=-1
418       dd=0
419    END IF
420
421    DO l = 1, llm
422       DO j=jj_begin-dd,jj_end+dd
423          DO i=ii_begin-dd,ii_end+dd
424             ij=(j-1)*iim+i
425             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
426             IF(comp) THEN
427                rhodz(ij,l) = m
428             ELSE
429                err = MAX(err,abs(m-rhodz(ij,l)))
430             END IF
431          ENDDO
432       ENDDO
433    ENDDO
434
435    IF(.NOT. comp) THEN
436       IF(err>1e-10) THEN
437          PRINT *, 'Discrepancy between ps and rhodz detected', err
438          STOP
439       ELSE
440!          PRINT *, 'No discrepancy between ps and rhodz detected'
441       END IF
442    END IF
443
444  END SUBROUTINE compute_rhodz
445
446  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
447    USE icosa
448    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
449    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
450    REAL(rstd), INTENT(IN) :: tau
451    LOGICAL, INTENT(INOUT) :: fluxt_zero
452    IF(fluxt_zero) THEN
453!       PRINT *, 'Accumulating fluxes (first)'
454       fluxt_zero=.FALSE.
455       hfluxt = tau*hflux
456       wfluxt = tau*wflux
457    ELSE
458!       PRINT *, 'Accumulating fluxes (next)'
459       hfluxt = hfluxt + tau*hflux
460       wfluxt = wfluxt + tau*wflux
461    END IF
462  END SUBROUTINE accumulate_fluxes
463 
464END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.