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

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 <euler>, <runge_kutta>, <leapfrog_matsuno>'
       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.

    CALL caldyn_BC(f_phis, f_wflux)

!$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