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

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

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