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

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

removed LeapFrog? & Adams-Bashforth from timeloop_gcm
File size: 16.8 KB
Line 
1MODULE timeloop_gcm_mod
2  USE transfert_mod
3  USE icosa
4  PRIVATE
5 
6  PUBLIC :: init_timeloop, timeloop
7
8  INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
9  INTEGER  :: itau_sync=10
10
11  TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0
12
13  TYPE(t_field),POINTER :: f_phis(:)
14  TYPE(t_field),POINTER :: f_q(:)
15  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
16  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
17  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
18  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
19  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
20  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
21  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
22  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:) 
23
24  INTEGER :: nb_stage, matsuno_period, scheme   
25
26  REAL(rstd),SAVE :: jD_cur, jH_cur
27  REAL(rstd),SAVE :: start_time
28
29CONTAINS
30 
31  SUBROUTINE init_timeloop
32  USE icosa
33  USE dissip_gcm_mod
34  USE caldyn_mod
35  USE etat0_mod
36  USE guided_mod
37  USE advect_tracer_mod
38  USE physics_mod
39  USE mpipara
40  USE omp_para
41  USE trace
42  USE transfert_mod
43  USE check_conserve_mod
44  USE ioipsl
45  IMPLICIT NONE
46
47    CHARACTER(len=255) :: scheme_name
48
49    CALL allocate_field(f_dps,field_t,type_real)
50    CALL allocate_field(f_du,field_u,type_real,llm)
51    CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
52! Model state at current time step (RK/MLF/Euler)
53    CALL allocate_field(f_phis,field_t,type_real)
54    CALL allocate_field(f_ps,field_t,type_real)
55    CALL allocate_field(f_u,field_u,type_real,llm)
56    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
57! Model state at previous time step (RK/MLF)
58    CALL allocate_field(f_psm1,field_t,type_real)
59    CALL allocate_field(f_um1,field_u,type_real,llm)
60    CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
61! Tracers
62    CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
63    CALL allocate_field(f_rhodz,field_t,type_real,llm)
64! Mass fluxes
65    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
66    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
67    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
68    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
69
70
71!----------------------------------------------------
72  IF (TRIM(time_style)=='lmd')  Then
73
74   day_step=180
75   CALL getin('day_step',day_step)
76
77   ndays=1
78   CALL getin('ndays',ndays)
79
80   dt = daysec/REAL(day_step)
81   itaumax = ndays*day_step
82
83   calend = 'earth_360d'
84   CALL getin('calend', calend)
85
86   day_ini = 0
87   CALL getin('day_ini',day_ini)
88
89   day_end = 0
90   CALL getin('day_end',day_end)
91
92   annee_ref = 1998
93   CALL getin('annee_ref',annee_ref)
94
95   start_time = 0
96   CALL getin('start_time',start_time) 
97
98   write_period=0
99   CALL getin('write_period',write_period)
100     
101   write_period=write_period/scale_factor
102   itau_out=FLOOR(write_period/dt)
103   
104   PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 
105
106  mois = 1 ; heure = 0.
107  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
108  jH_ref = jD_ref - int(jD_ref) 
109  jD_ref = int(jD_ref) 
110 
111  CALL ioconf_startdate(INT(jD_ref),jH_ref) 
112  write(*,*)'annee_ref, mois, day_ref, heure, jD_ref'
113  write(*,*)annee_ref, mois, day_ref, heure, jD_ref
114  write(*,*)"ndays,day_step,itaumax,dt======>"
115  write(*,*)ndays,day_step,itaumax,dt
116  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
117  write(*,*)'jD_ref+jH_ref,an, mois, jour, heure'
118  write(*,*)jD_ref+jH_ref,an, mois, jour, heure
119  day_end = day_ini + ndays 
120 END IF 
121!----------------------------------------------------
122
123
124
125
126
127    scheme_name='runge_kutta'
128    CALL getin('scheme',scheme_name)
129
130    SELECT CASE (TRIM(scheme_name))
131      CASE('euler')
132        scheme=euler
133        nb_stage=1
134      CASE ('runge_kutta')
135        scheme=rk4
136        nb_stage=4
137      CASE ('leapfrog_matsuno')
138        scheme=mlf
139        matsuno_period=5
140        CALL getin('matsuno_period',matsuno_period)
141        nb_stage=matsuno_period+1
142     ! Model state 2 time steps ago (MLF)
143        CALL allocate_field(f_psm2,field_t,type_real)
144        CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
145        CALL allocate_field(f_um2,field_u,type_real,llm)
146    CASE default
147       PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
148            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
149       STOP
150    END SELECT
151
152
153    CALL init_dissip
154    CALL init_caldyn
155    CALL init_guided
156    CALL init_advect_tracer
157    CALL init_check_conserve
158    CALL init_physics
159    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
160
161    CALL transfert_request(f_phis,req_i0) 
162    CALL transfert_request(f_phis,req_i1) 
163    CALL writefield("phis",f_phis,once=.TRUE.)
164
165    CALL init_message(f_ps,req_i0,req_ps0)
166    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
167    CALL init_message(f_u,req_e0_vect,req_u0)
168    CALL init_message(f_q,req_i0,req_q0)
169
170  END SUBROUTINE init_timeloop
171 
172  SUBROUTINE timeloop
173  USE icosa
174  USE dissip_gcm_mod
175  USE caldyn_mod
176  USE etat0_mod
177  USE guided_mod
178  USE advect_tracer_mod
179  USE physics_mod
180  USE mpipara
181  USE omp_para
182  USE trace
183  USE transfert_mod
184  USE check_conserve_mod
185  IMPLICIT NONE 
186    REAL(rstd),POINTER :: phis(:)
187    REAL(rstd),POINTER :: q(:,:,:)
188    REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
189    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
190    REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
191    REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
192    REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
193    REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
194    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
195
196    INTEGER :: ind
197    INTEGER :: it,i,j,n,  stage
198    CHARACTER(len=255) :: scheme_name
199    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
200    LOGICAL, PARAMETER :: check=.FALSE.
201
202    CALL caldyn_BC(f_phis, f_wflux)
203
204!$OMP BARRIER
205  DO ind=1,ndomain
206     CALL swap_dimensions(ind)
207     CALL swap_geometry(ind)
208     rhodz=f_rhodz(ind); ps=f_ps(ind)
209     CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
210  END DO
211  fluxt_zero=.TRUE.
212
213  DO it=0,itaumax
214    IF (MOD(it,itau_sync)==0) THEN
215      CALL send_message(f_ps,req_ps0)
216      CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
217      CALL send_message(f_u,req_u0)
218      CALL send_message(f_q,req_q0) 
219      CALL wait_message(req_ps0)
220      CALL wait_message(req_theta_rhodz0) 
221      CALL wait_message(req_u0)
222      CALL wait_message(req_q0) 
223    ENDIF
224   
225!    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
226    IF (mod(it,itau_out)==0 ) THEN
227      CALL writefield("q",f_q)
228      CALL update_time_counter(dt*it)
229      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
230    ENDIF
231   
232    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
233
234    DO stage=1,nb_stage
235       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
236            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
237            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
238       SELECT CASE (scheme)
239       CASE(euler)
240          CALL euler_scheme(.TRUE.)
241       CASE (rk4)
242          CALL rk_scheme(stage)
243       CASE (mlf)
244          CALL  leapfrog_matsuno_scheme(stage)
245       CASE DEFAULT
246          STOP
247       END SELECT
248    END DO
249
250    IF (MOD(it+1,itau_dissip)==0) THEN
251      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
252      CALL euler_scheme(.FALSE.)
253    ENDIF
254
255    IF(MOD(it+1,itau_adv)==0) THEN
256
257       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
258       fluxt_zero=.TRUE.
259
260       ! FIXME : check that rhodz is consistent with ps
261       IF (check) THEN
262         DO ind=1,ndomain
263            CALL swap_dimensions(ind)
264            CALL swap_geometry(ind)
265            rhodz=f_rhodz(ind); ps=f_ps(ind);
266            CALL compute_rhodz(.FALSE., ps, rhodz)   
267         END DO
268       ENDIF
269
270    END IF
271
272
273
274!----------------------------------------------------
275!    jD_cur = jD_ref + day_ini - day_ref + it/day_step
276!    jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
277!    jD_cur = jD_cur + int(jH_cur)
278!    jH_cur = jH_cur - int(jH_cur)
279!    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
280
281!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
282    ENDDO
283
284 
285  CONTAINS
286
287    SUBROUTINE Euler_scheme(with_dps)
288    IMPLICIT NONE
289    LOGICAL :: with_dps
290    INTEGER :: ind
291    INTEGER :: i,j,ij,l
292    CALL trace_start("Euler_scheme") 
293
294    DO ind=1,ndomain
295       CALL swap_dimensions(ind)
296       CALL swap_geometry(ind)
297       IF(with_dps) THEN
298         ps=f_ps(ind) ; dps=f_dps(ind) ; 
299
300         IF (omp_first) THEN
301           DO j=jj_begin,jj_end
302             DO i=ii_begin,ii_end
303               ij=(j-1)*iim+i
304               ps(ij)=ps(ij)+dt*dps(ij)
305             ENDDO
306           ENDDO
307         ENDIF
308         
309         hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
310         wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
311         CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
312       END IF
313       
314       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
315       du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
316
317       DO l=ll_begin,ll_end
318         DO j=jj_begin,jj_end
319           DO i=ii_begin,ii_end
320             ij=(j-1)*iim+i
321             u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
322             u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
323             u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
324             theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
325           ENDDO
326         ENDDO
327       ENDDO
328    ENDDO
329
330    CALL trace_end("Euler_scheme") 
331
332    END SUBROUTINE Euler_scheme
333
334    SUBROUTINE RK_scheme(stage)
335      USE caldyn_gcm_mod
336      IMPLICIT NONE
337      INTEGER :: ind, stage
338      REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
339      REAL(rstd) :: tau
340      INTEGER :: i,j,ij,l
341 
342      CALL trace_start("RK_scheme") 
343
344      tau = dt*coef(stage)
345
346      DO ind=1,ndomain
347         CALL swap_dimensions(ind)
348         CALL swap_geometry(ind)
349         ps=f_ps(ind)   
350         psm1=f_psm1(ind) 
351         dps=f_dps(ind) 
352         
353         IF (stage==1) THEN ! first stage : save model state in XXm1
354           IF (omp_first) THEN
355             DO j=jj_begin,jj_end
356               DO i=ii_begin,ii_end
357                 ij=(j-1)*iim+i
358                 psm1(ij)=ps(ij)
359               ENDDO
360             ENDDO
361           ENDIF
362         END IF
363         
364         ! updates are of the form x1 := x0 + tau*f(x1)
365           IF (omp_first) THEN
366             DO j=jj_begin,jj_end
367               DO i=ii_begin,ii_end
368                 ij=(j-1)*iim+i
369                 ps(ij)=psm1(ij)+tau*dps(ij)
370               ENDDO
371             ENDDO
372           ENDIF
373       
374       ENDDO
375
376       CALL send_message(f_ps,req_ps)     
377
378
379     
380      DO ind=1,ndomain
381         CALL swap_dimensions(ind)
382         CALL swap_geometry(ind)
383         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
384         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
385         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
386         
387         IF (stage==1) THEN ! first stage : save model state in XXm1
388             
389         DO l=ll_begin,ll_end
390             DO j=jj_begin,jj_end
391               DO i=ii_begin,ii_end
392                 ij=(j-1)*iim+i
393                 um1(ij+u_right,l)=u(ij+u_right,l)
394                 um1(ij+u_lup,l)=u(ij+u_lup,l)
395                 um1(ij+u_ldown,l)=u(ij+u_ldown,l)
396                 theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
397               ENDDO
398             ENDDO
399           ENDDO
400
401         END IF
402         ! updates are of the form x1 := x0 + tau*f(x1)
403         
404         DO l=ll_begin,ll_end
405           DO j=jj_begin,jj_end
406             DO i=ii_begin,ii_end
407               ij=(j-1)*iim+i
408               u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
409               u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
410               u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
411               theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
412             ENDDO
413           ENDDO
414         ENDDO
415
416         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
417            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
418            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
419            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
420         END IF
421      END DO
422     
423      CALL trace_end("RK_scheme")
424     
425    END SUBROUTINE RK_scheme
426
427    SUBROUTINE leapfrog_matsuno_scheme(stage)
428    IMPLICIT NONE
429    INTEGER :: ind, stage
430    REAL :: tau
431
432      CALL trace_start("leapfrog_matsuno_scheme")
433   
434      tau = dt/nb_stage
435      DO ind=1,ndomain
436        CALL swap_dimensions(ind)
437        CALL swap_geometry(ind)
438
439        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
440        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
441        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
442        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
443
444       
445        IF (stage==1) THEN
446          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
447          ps(:)=psm1(:)+tau*dps(:)
448          u(:,:)=um1(:,:)+tau*du(:,:)
449          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
450
451        ELSE IF (stage==2) THEN
452
453          ps(:)=psm1(:)+tau*dps(:)
454          u(:,:)=um1(:,:)+tau*du(:,:)
455          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
456
457          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
458          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
459
460        ELSE
461
462          ps(:)=psm2(:)+2*tau*dps(:)
463          u(:,:)=um2(:,:)+2*tau*du(:,:)
464          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
465
466          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
467          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
468
469        ENDIF
470     
471      ENDDO
472      CALL trace_end("leapfrog_matsuno_scheme")
473     
474    END SUBROUTINE leapfrog_matsuno_scheme 
475         
476  END SUBROUTINE timeloop   
477
478  SUBROUTINE compute_rhodz(comp, ps, rhodz)
479    USE icosa
480    USE disvert_mod
481    USE omp_para
482    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
483    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
484    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
485    REAL(rstd) :: m, err
486    INTEGER :: l,i,j,ij,dd
487    err=0.
488
489    dd=0
490
491    DO l = ll_begin, ll_end
492       DO j=jj_begin-dd,jj_end+dd
493          DO i=ii_begin-dd,ii_end+dd
494             ij=(j-1)*iim+i
495             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
496             IF(comp) THEN
497                rhodz(ij,l) = m
498             ELSE
499                err = MAX(err,abs(m-rhodz(ij,l)))
500             END IF
501          ENDDO
502       ENDDO
503    ENDDO
504
505    IF(.NOT. comp) THEN
506       IF(err>1e-10) THEN
507          PRINT *, 'Discrepancy between ps and rhodz detected', err
508          STOP
509       ELSE
510!          PRINT *, 'No discrepancy between ps and rhodz detected'
511       END IF
512    END IF
513
514  END SUBROUTINE compute_rhodz
515
516  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
517    USE icosa
518    USE omp_para
519    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
520    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
521    REAL(rstd), INTENT(IN) :: tau
522    LOGICAL, INTENT(INOUT) :: fluxt_zero
523    INTEGER :: l,i,j,ij
524
525    IF(fluxt_zero) THEN
526
527       fluxt_zero=.FALSE.
528
529       DO l=ll_begin,ll_end
530         DO j=jj_begin-1,jj_end+1
531           DO i=ii_begin-1,ii_end+1
532             ij=(j-1)*iim+i
533             hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
534             hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
535             hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
536           ENDDO
537         ENDDO
538       ENDDO
539
540       DO l=ll_begin,ll_endp1
541         DO j=jj_begin,jj_end
542           DO i=ii_begin,ii_end
543             ij=(j-1)*iim+i
544             wfluxt(ij,l) = tau*wflux(ij,l)
545           ENDDO
546         ENDDO
547       ENDDO
548    ELSE
549
550       DO l=ll_begin,ll_end
551         DO j=jj_begin-1,jj_end+1
552           DO i=ii_begin-1,ii_end+1
553             ij=(j-1)*iim+i
554             hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
555             hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
556             hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
557           ENDDO
558         ENDDO
559       ENDDO
560
561       DO l=ll_begin,ll_endp1
562         DO j=jj_begin,jj_end
563           DO i=ii_begin,ii_end
564             ij=(j-1)*iim+i
565             wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
566           ENDDO
567         ENDDO
568       ENDDO
569
570    END IF
571
572  END SUBROUTINE accumulate_fluxes
573 
574END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.