MODULE timeloop_gcm_mod USE transfert_mod USE icosa PRIVATE PUBLIC :: init_timeloop, timeloop INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3 INTEGER :: itau_sync=10 TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0 TYPE(t_field),POINTER :: f_phis(:) TYPE(t_field),POINTER :: f_q(:) TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:) TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:) TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:) TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:) TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:) TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:) TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:) TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:) INTEGER :: nb_stage, matsuno_period, scheme REAL(rstd),SAVE :: jD_cur, jH_cur REAL(rstd),SAVE :: start_time CONTAINS SUBROUTINE init_timeloop USE icosa USE dissip_gcm_mod USE caldyn_mod USE etat0_mod USE guided_mod USE advect_tracer_mod USE physics_mod USE mpipara USE omp_para USE trace USE transfert_mod USE check_conserve_mod USE ioipsl IMPLICIT NONE CHARACTER(len=255) :: scheme_name CALL allocate_field(f_dps,field_t,type_real) CALL allocate_field(f_du,field_u,type_real,llm) CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm) ! Model state at current time step (RK/MLF/Euler) CALL allocate_field(f_phis,field_t,type_real) CALL allocate_field(f_ps,field_t,type_real) CALL allocate_field(f_u,field_u,type_real,llm) CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) ! Model state at previous time step (RK/MLF) CALL allocate_field(f_psm1,field_t,type_real) CALL allocate_field(f_um1,field_u,type_real,llm) CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) ! Tracers CALL allocate_field(f_q,field_t,type_real,llm,nqtot) CALL allocate_field(f_rhodz,field_t,type_real,llm) ! Mass fluxes CALL allocate_field(f_hflux,field_u,type_real,llm) ! instantaneous mass fluxes CALL allocate_field(f_wflux,field_t,type_real,llm+1) ! computed by caldyn CALL allocate_field(f_hfluxt,field_u,type_real,llm) ! mass "fluxes" accumulated in time CALL allocate_field(f_wfluxt,field_t,type_real,llm+1) !---------------------------------------------------- IF (TRIM(time_style)=='lmd') Then day_step=180 CALL getin('day_step',day_step) ndays=1 CALL getin('ndays',ndays) dt = daysec/REAL(day_step) itaumax = ndays*day_step calend = 'earth_360d' CALL getin('calend', calend) day_ini = 0 CALL getin('day_ini',day_ini) day_end = 0 CALL getin('day_end',day_end) annee_ref = 1998 CALL getin('annee_ref',annee_ref) start_time = 0 CALL getin('start_time',start_time) write_period=0 CALL getin('write_period',write_period) write_period=write_period/scale_factor itau_out=FLOOR(write_period/dt) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out mois = 1 ; heure = 0. call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref) jH_ref = jD_ref - int(jD_ref) jD_ref = int(jD_ref) CALL ioconf_startdate(INT(jD_ref),jH_ref) write(*,*)'annee_ref, mois, day_ref, heure, jD_ref' write(*,*)annee_ref, mois, day_ref, heure, jD_ref write(*,*)"ndays,day_step,itaumax,dt======>" write(*,*)ndays,day_step,itaumax,dt call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure) write(*,*)'jD_ref+jH_ref,an, mois, jour, heure' write(*,*)jD_ref+jH_ref,an, mois, jour, heure day_end = day_ini + ndays END IF !---------------------------------------------------- scheme_name='runge_kutta' CALL getin('scheme',scheme_name) SELECT CASE (TRIM(scheme_name)) CASE('euler') scheme=euler nb_stage=1 CASE ('runge_kutta') scheme=rk4 nb_stage=4 CASE ('leapfrog_matsuno') scheme=mlf matsuno_period=5 CALL getin('matsuno_period',matsuno_period) nb_stage=matsuno_period+1 ! Model state 2 time steps ago (MLF) CALL allocate_field(f_psm2,field_t,type_real) CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) CALL allocate_field(f_um2,field_u,type_real,llm) CASE default PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name), & ' > options are , , ' STOP END SELECT CALL init_dissip CALL init_caldyn CALL init_guided CALL init_advect_tracer CALL init_check_conserve CALL init_physics CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CALL transfert_request(f_phis,req_i0) CALL transfert_request(f_phis,req_i1) CALL writefield("phis",f_phis,once=.TRUE.) CALL init_message(f_ps,req_i0,req_ps0) CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0) CALL init_message(f_u,req_e0_vect,req_u0) CALL init_message(f_q,req_i0,req_q0) END SUBROUTINE init_timeloop SUBROUTINE timeloop USE icosa USE dissip_gcm_mod USE caldyn_mod USE etat0_mod USE guided_mod USE advect_tracer_mod USE physics_mod USE mpipara USE omp_para USE trace USE transfert_mod USE check_conserve_mod IMPLICIT NONE REAL(rstd),POINTER :: phis(:) REAL(rstd),POINTER :: q(:,:,:) REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:) REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:) REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:) REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:) REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:) REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:) REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) INTEGER :: ind INTEGER :: it,i,j,n, stage CHARACTER(len=255) :: scheme_name LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time LOGICAL, PARAMETER :: check=.FALSE. !$OMP BARRIER DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) rhodz=f_rhodz(ind); ps=f_ps(ind) CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps END DO fluxt_zero=.TRUE. DO it=0,itaumax IF (MOD(it,itau_sync)==0) THEN CALL send_message(f_ps,req_ps0) CALL send_message(f_theta_rhodz,req_theta_rhodz0) CALL send_message(f_u,req_u0) CALL send_message(f_q,req_q0) CALL wait_message(req_ps0) CALL wait_message(req_theta_rhodz0) CALL wait_message(req_u0) CALL wait_message(req_q0) ENDIF ! IF (is_mpi_root) PRINT *,"It No :",It," t :",dt*It IF (mod(it,itau_out)==0 ) THEN CALL writefield("q",f_q) CALL update_time_counter(dt*it) CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) ENDIF CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) DO stage=1,nb_stage CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & f_phis,f_ps,f_theta_rhodz,f_u, f_q, & f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) SELECT CASE (scheme) CASE(euler) CALL euler_scheme(.TRUE.) CASE (rk4) CALL rk_scheme(stage) CASE (mlf) CALL leapfrog_matsuno_scheme(stage) ! CASE ('leapfrog') ! CALL leapfrog_scheme ! ! CASE ('adam_bashforth') ! CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) ! CALL adam_bashforth_scheme CASE DEFAULT STOP END SELECT END DO IF (MOD(it+1,itau_dissip)==0) THEN CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) CALL euler_scheme(.FALSE.) ENDIF IF(MOD(it+1,itau_adv)==0) THEN CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step fluxt_zero=.TRUE. ! FIXME : check that rhodz is consistent with ps IF (check) THEN DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) rhodz=f_rhodz(ind); ps=f_ps(ind); CALL compute_rhodz(.FALSE., ps, rhodz) END DO ENDIF END IF !---------------------------------------------------- ! jD_cur = jD_ref + day_ini - day_ref + it/day_step ! jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step) ! jD_cur = jD_cur + int(jH_cur) ! jH_cur = jH_cur - int(jH_cur) ! CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q) ! CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q) ENDDO CONTAINS SUBROUTINE Euler_scheme(with_dps) IMPLICIT NONE LOGICAL :: with_dps INTEGER :: ind INTEGER :: i,j,ij,l CALL trace_start("Euler_scheme") DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) IF(with_dps) THEN ps=f_ps(ind) ; dps=f_dps(ind) ; IF (omp_first) THEN DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i ps(ij)=ps(ij)+dt*dps(ij) ENDDO ENDDO ENDIF hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind)) END IF u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l) u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l) u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l) theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l) ENDDO ENDDO ENDDO ENDDO CALL trace_end("Euler_scheme") END SUBROUTINE Euler_scheme SUBROUTINE RK_scheme(stage) USE caldyn_gcm_mod IMPLICIT NONE INTEGER :: ind, stage REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /) REAL(rstd) :: tau INTEGER :: i,j,ij,l CALL trace_start("RK_scheme") tau = dt*coef(stage) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) psm1=f_psm1(ind) dps=f_dps(ind) IF (stage==1) THEN ! first stage : save model state in XXm1 IF (omp_first) THEN DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i psm1(ij)=ps(ij) ENDDO ENDDO ENDIF END IF ! updates are of the form x1 := x0 + tau*f(x1) IF (omp_first) THEN DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i ps(ij)=psm1(ij)+tau*dps(ij) ENDDO ENDDO ENDIF ENDDO CALL send_message(f_ps,req_ps) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) IF (stage==1) THEN ! first stage : save model state in XXm1 DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i um1(ij+u_right,l)=u(ij+u_right,l) um1(ij+u_lup,l)=u(ij+u_lup,l) um1(ij+u_ldown,l)=u(ij+u_ldown,l) theta_rhodzm1(ij,l)=theta_rhodz(ij,l) ENDDO ENDDO ENDDO END IF ! updates are of the form x1 := x0 + tau*f(x1) DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l) u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l) u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l) theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l) ENDDO ENDDO ENDDO IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind)) END IF END DO CALL trace_end("RK_scheme") END SUBROUTINE RK_scheme SUBROUTINE leapfrog_scheme IMPLICIT NONE INTEGER :: ind CALL trace_start("leapfrog_scheme") DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind) dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) IF (it==0) THEN psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) ps(:)=ps(:)+dt*dps(:) u(:,:)=u(:,:)+dt*du(:,:) theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:) psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) ELSE ps(:)=psm2(:)+2*dt*dps(:) u(:,:)=um2(:,:)+2*dt*du(:,:) theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:) psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) ENDIF ENDDO CALL trace_end("leapfrog_scheme") END SUBROUTINE leapfrog_scheme SUBROUTINE leapfrog_matsuno_scheme(stage) IMPLICIT NONE INTEGER :: ind, stage REAL :: tau CALL trace_start("leapfrog_matsuno_scheme") tau = dt/nb_stage DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind) dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) IF (stage==1) THEN psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ps(:)=psm1(:)+tau*dps(:) u(:,:)=um1(:,:)+tau*du(:,:) theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) ELSE IF (stage==2) THEN ps(:)=psm1(:)+tau*dps(:) u(:,:)=um1(:,:)+tau*du(:,:) theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) ELSE ps(:)=psm2(:)+2*tau*dps(:) u(:,:)=um2(:,:)+2*tau*du(:,:) theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:) psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) ENDIF ENDDO CALL trace_end("leapfrog_matsuno_scheme") END SUBROUTINE leapfrog_matsuno_scheme SUBROUTINE adam_bashforth_scheme IMPLICIT NONE INTEGER :: ind CALL trace_start("adam_bashforth_scheme") DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind) dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind) IF (it==0) THEN dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:) dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:) ENDIF ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12 u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12 theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12 dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:) dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:) ENDDO CALL trace_end("adam_bashforth_scheme") END SUBROUTINE adam_bashforth_scheme END SUBROUTINE timeloop SUBROUTINE compute_rhodz(comp, ps, rhodz) USE icosa USE disvert_mod USE omp_para LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check REAL(rstd), INTENT(IN) :: ps(iim*jjm) REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm) REAL(rstd) :: m, err INTEGER :: l,i,j,ij,dd err=0. dd=0 DO l = ll_begin, ll_end DO j=jj_begin-dd,jj_end+dd DO i=ii_begin-dd,ii_end+dd ij=(j-1)*iim+i m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g IF(comp) THEN rhodz(ij,l) = m ELSE err = MAX(err,abs(m-rhodz(ij,l))) END IF ENDDO ENDDO ENDDO IF(.NOT. comp) THEN IF(err>1e-10) THEN PRINT *, 'Discrepancy between ps and rhodz detected', err STOP ELSE ! PRINT *, 'No discrepancy between ps and rhodz detected' END IF END IF END SUBROUTINE compute_rhodz SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) USE icosa USE omp_para REAL(rstd), INTENT(IN) :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1) REAL(rstd), INTENT(IN) :: tau LOGICAL, INTENT(INOUT) :: fluxt_zero INTEGER :: l,i,j,ij IF(fluxt_zero) THEN fluxt_zero=.FALSE. DO l=ll_begin,ll_end DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l) hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l) hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l) ENDDO ENDDO ENDDO DO l=ll_begin,ll_endp1 DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i wfluxt(ij,l) = tau*wflux(ij,l) ENDDO ENDDO ENDDO ELSE DO l=ll_begin,ll_end DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l) hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l) hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l) ENDDO ENDDO ENDDO DO l=ll_begin,ll_endp1 DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) ENDDO ENDDO ENDDO END IF END SUBROUTINE accumulate_fluxes END MODULE timeloop_gcm_mod