MODULE caldyn_gcm_mod
  USE icosa

  TYPE(t_field),POINTER :: f_out_u(:)
  REAL(rstd),POINTER :: out_u(:,:)

  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)

  PUBLIC init_caldyn, caldyn, write_output_fields

  INTEGER :: caldyn_hydrostat

CONTAINS
  
  SUBROUTINE init_caldyn
    USE icosa
    USE exner_mod
    IMPLICIT NONE
    CHARACTER(len=255) :: def
  
    def='direct'
    CALL getin('caldyn_exner',def)
    SELECT CASE(TRIM(def))
    CASE('lmdz')
       caldyn_exner=1
    CASE('direct')
       caldyn_exner=2
    CASE DEFAULT
       PRINT*,'Bad selector for variable caldyn_exner : <', TRIM(def),'> options are <lmdz>, <direct>'
       STOP
    END SELECT

    def='direct'
    CALL getin('caldyn_hydrostat',def)
    SELECT CASE(TRIM(def))
    CASE('lmdz')
       caldyn_hydrostat=1
    CASE('direct')
       caldyn_hydrostat=2
    CASE DEFAULT
       PRINT*,'Bad selector for variable caldyn_hydrostat : <', TRIM(def),'> options are <lmdz>, <direct>'
       STOP
    END SELECT
    
    CALL allocate_caldyn
  
  END SUBROUTINE init_caldyn
 
  SUBROUTINE allocate_caldyn
  USE icosa
  IMPLICIT NONE

  CALL allocate_field(f_out_u,field_u,type_real,llm) 
  
  CALL allocate_field(f_buf_i,field_t,type_real,llm)
  CALL allocate_field(f_buf_p,field_t,type_real,llm+1)
  CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
  CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
  CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
  CALL allocate_field(f_buf_v,field_z,type_real,llm)
  CALL allocate_field(f_buf_s,field_t,type_real)

  END SUBROUTINE allocate_caldyn
   
  SUBROUTINE check_mass_conservation(f_ps,f_dps)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_dps(:)
    REAL(rstd),POINTER :: ps(:)
    REAL(rstd),POINTER :: dps(:)
    REAL(rstd) :: mass_tot,dmass_tot
    INTEGER :: ind,i,j,ij
    
    mass_tot=0
    dmass_tot=0
    
    CALL transfert_request(f_dps,req_i1)
    CALL transfert_request(f_ps,req_i1)

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)

      ps=f_ps(ind)
      dps=f_dps(ind)

      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          IF (domain(ind)%own(i,j)) THEN
            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
          ENDIF
        ENDDO
      ENDDO
   
    ENDDO
    PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot        

  END SUBROUTINE check_mass_conservation
  
  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
  USE icosa
  USE vorticity_mod
  USE kinetic_mod
  USE theta2theta_rhodz_mod
  IMPLICIT NONE
  INTEGER,INTENT(IN)    :: it
  TYPE(t_field),POINTER :: f_phis(:)
  TYPE(t_field),POINTER :: f_ps(:)
  TYPE(t_field),POINTER :: f_theta_rhodz(:)
  TYPE(t_field),POINTER :: f_u(:)
  TYPE(t_field),POINTER :: f_q(:)
  TYPE(t_field),POINTER :: f_dps(:)
  TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
  TYPE(t_field),POINTER :: f_du(:)

  REAL(rstd),POINTER :: phis(:)
  REAL(rstd),POINTER :: ps(:)
  REAL(rstd),POINTER :: theta_rhodz(:,:)
  REAL(rstd),POINTER :: u(:,:)
  REAL(rstd),POINTER :: dps(:)
  REAL(rstd),POINTER :: dtheta_rhodz(:,:)
  REAL(rstd),POINTER :: du(:,:)
  INTEGER :: ind,ij

  
    CALL transfert_request(f_phis,req_i1) 
    CALL transfert_request(f_ps,req_i1) 
    CALL transfert_request(f_theta_rhodz,req_i1) 
    CALL transfert_request(f_u,req_e1)
    
    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      
      out_u=f_out_u(ind) 
      phis=f_phis(ind)
      ps=f_ps(ind)
      theta_rhodz=f_theta_rhodz(ind)
      u=f_u(ind)
      dps=f_dps(ind)
      dtheta_rhodz=f_dtheta_rhodz(ind)
      du=f_du(ind)
!$OMP PARALLEL DEFAULT(SHARED)
      CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
!$OMP END PARALLEL
    ENDDO

    IF (mod(it,itau_out)==0 ) THEN
       PRINT *,'CALL write_output_fields'
       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
    END IF

!    CALL check_mass_conservation(f_ps,f_dps)

  END SUBROUTINE caldyn
  
  
  SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
  USE icosa
  USE disvert_mod
  USE exner_mod
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
    REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm)
    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
    REAL(rstd),INTENT(OUT):: dtheta_rhodz(iim*jjm,llm)
    REAL(rstd),INTENT(OUT):: dps(iim*jjm)
    
    INTEGER :: i,j,ij,l
    REAL(rstd) :: ww,uu 
    REAL(rstd) :: delta
    REAL(rstd) :: etav,hv, du2

    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)  ! pression
    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:)   ! Exner function
    REAL(rstd),ALLOCATABLE,SAVE :: pks(:)
 ! Intermediate variable to compute exner function
    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:)
    REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:)
 !   
    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
    REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:)   ! mass density   
    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux   
    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux   
    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)    ! mass flux convergence
    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)       ! vertical velocity      
    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)       ! potential velocity  
    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)   ! bernouilli term
    
    LOGICAL,SAVE :: first=.TRUE.
!$OMP THREADPRIVATE(first)
        
!$OMP BARRIER      
!$OMP MASTER  
!    IF (first) THEN
      ALLOCATE(theta(iim*jjm,llm))  ! potential temperature
      ALLOCATE(p(iim*jjm,llm+1))  ! pression
      ALLOCATE(pk(iim*jjm,llm))   ! Exner function
      ALLOCATE(pks(iim*jjm))
      ALLOCATE(alpha(iim*jjm,llm))
      ALLOCATE(beta(iim*jjm,llm))
      ALLOCATE(phi(iim*jjm,llm))    ! geopotential
      ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
      ALLOCATE(Fe(3*iim*jjm,llm))   ! mass flux   
      ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
      ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
      ALLOCATE(w(iim*jjm,llm))       ! vertical velocity      
      ALLOCATE(qv(2*iim*jjm,llm))       ! potential velocity  
      ALLOCATE(berni(iim*jjm,llm))   ! bernouilli term
!      first=.FALSE.
!    ENDIF
!$OMP END MASTER
!$OMP BARRIER        
 !    du(:,:)=0
!    theta=1e10
!    p=1e10
!    pk=1e10
!    pks=1e10
!    alpha=1e10
!    beta=1e10
!    phi=1e10
!    mass=1e10
!    rhodz=1e10
!    Fe=1e10
!    Ftheta=1e10
!    convm=1e10
!    w=1e10
!    qv=1e10
!    berni=1e10
    
 !!! Compute pressure

!    PRINT *, 'Computing pressure'
    DO    l    = 1, llm+1
!$OMP DO
      DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
          ij=(j-1)*iim+i
          p(ij,l) = ap(l) + bp(l) * ps(ij)
        ENDDO
      ENDDO
    ENDDO

!!! Compute mass
   ! PRINT *, 'Computing mass, theta'
   DO l = 1, llm
!$OMP DO
      DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
            ij=(j-1)*iim+i
            rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
         ENDDO
      ENDDO
   ENDDO

!!! Compute shallow-water potential vorticity
  DO l = 1,llm
!$OMP DO
    DO j=jj_begin-1,jj_end+1
      DO i=ii_begin-1,ii_end+1
        ij=(j-1)*iim+i
        
          etav= 1./Av(ij+z_up)*(  ne(ij,rup)        * u(ij+u_rup,l)        * de(ij+u_rup)         &
                                + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
                                - ne(ij,lup)        * u(ij+u_lup,l)        * de(ij+u_lup) )                               

          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
     
          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
          
          etav = 1./Av(ij+z_down)*(  ne(ij,ldown)         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
                                   + ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
                                   - ne(ij,rdown)         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
       
          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
                       
          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv

        ENDDO
      ENDDO
    ENDDO


  DO l = 1, llm
!!!  Compute mass and theta fluxes
!$OMP DO
    DO j=jj_begin-1,jj_end+1
      DO i=ii_begin-1,ii_end+1
        ij=(j-1)*iim+i
        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
        Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)

        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l)
        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l)
        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l)
      ENDDO
    ENDDO

!!! compute horizontal divergence of fluxes
!$OMP DO
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i  
        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
                                ne(ij,rup)*Fe(ij+u_rup,l)       +  &  
                                ne(ij,lup)*Fe(ij+u_lup,l)       +  &  
                                ne(ij,left)*Fe(ij+u_left,l)     +  &  
                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &  
                                ne(ij,rdown)*Fe(ij+u_rdown,l))    

        ! signe ? attention d (rho theta dz)
        ! dtheta_rhodz =  -div(flux.theta)
        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &  
                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  &  
                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  &  
                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  &  
                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
      ENDDO
    ENDDO
  ENDDO
   
!!! cumulate mass flux convergence from top to bottom
  DO  l = llm-1, 1, -1
!$OMP DO
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
      ENDDO
    ENDDO
  ENDDO
  
!!! Compute vertical mass flux
  DO l = 1,llm-1
!$OMP DO
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
        ! => w>0 for upward transport
        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
      ENDDO
    ENDDO
  ENDDO

! compute dps, vertical mass flux at the surface = 0
!$OMP DO
  DO j=jj_begin,jj_end
    DO i=ii_begin,ii_end
      ij=(j-1)*iim+i
      w(ij,1)  = 0.
      ! dps/dt = -int(div flux)dz
      dps(ij)=-convm(ij,1) * g 
    ENDDO
  ENDDO

!!! Compute potential vorticity (Coriolis) contribution to du
  DO l=1,llm
!$OMP DO
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i

       du(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))/de(ij+u_right) *      &  
                      ( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+                           &
                        wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+                           &
                        wee(ij+u_right,3,1)*Fe(ij+u_left,l)+                          &
                        wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+                         &
                        wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+                         & 
                        wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+                 &
                        wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+                 &
                        wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+                 &
                        wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+                   &
                        wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) )


       du(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))/de(ij+u_lup) *         &  
                      ( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+                        &
                        wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+                       &
                        wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+                       &
                        wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+                       &
                        wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+                         & 
                        wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+                 &
                        wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+                   &
                        wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+                   &
                        wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+                  &
                        wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) )


       du(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))/de(ij+u_ldown) *   &  
                      ( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+                      &
                        wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+                      &
                        wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+                        &
                        wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+                        &
                        wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+                       & 
                        wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+                &
                        wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+               &
                        wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+              &
                        wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+              &
                        wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) )

     
      ENDDO
    ENDDO
  ENDDO
  
!!! Compute Exner function
!    PRINT *, 'Computing Exner'
    CALL compute_exner(ps,p,pks,pk,1) 

!!! Compute geopotential

  ! for first layer
!$OMP DO
   DO j=jj_begin-1,jj_end+1
     DO i=ii_begin-1,ii_end+1
       ij=(j-1)*iim+i
       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
     ENDDO
   ENDDO
          
  ! for other layers
   DO l = 2, llm
!$OMP DO
     DO j=jj_begin-1,jj_end+1
       DO i=ii_begin-1,ii_end+1
         ij=(j-1)*iim+i
         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
                                       * (  pk(ij,l-1) -  pk(ij,l)    )
       ENDDO
     ENDDO
   ENDDO

   
!!! Compute bernouilli term = Kinetic Energy + geopotential
  DO l=1,llm
!$OMP DO
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i

        berni(ij,l) =  phi(ij,l) &                                                          
                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
       
      ENDDO
    ENDDO
  ENDDO
  
 
!!! gradients of Bernoulli and Exner functions
  DO l=1,llm
!$OMP DO
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        
        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
                             *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
        
        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l)
        
        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
        
        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l)
                
        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
        
        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l)

      ENDDO
    ENDDO
  ENDDO
 
!!! contributions of vertical advection to du, dtheta
 
  DO l=1,llm-1
!$OMP DO
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
         ! ww>0 <=> upward transport

        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww

        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1))
        uu  = u(ij+u_right,l+1) - u(ij+u_right,l)  
        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
        du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))

        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_lup,l+1))
        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l)  
        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
        du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))

        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_ldown,l+1))
        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l)  
        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
        du(ij+u_ldown, l+1 ) = du(ij+u_ldown,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_ldown,l+1)))

      ENDDO
    ENDDO
  ENDDO      
   
!!$OMP BARRIER
!!$OMP MASTER  
    DEALLOCATE(theta)  ! potential temperature
    DEALLOCATE(p)  ! pression
    DEALLOCATE(pk)   ! Exner function
    DEALLOCATE(pks)
    DEALLOCATE(alpha)
    DEALLOCATE(beta)
    DEALLOCATE(phi)    ! geopotential
!    DEALLOCATE(mass)   ! mass    
    DEALLOCATE(rhodz)   ! mass density   
    DEALLOCATE(Fe)   ! mass flux   
    DEALLOCATE(Ftheta) ! theta flux   
    DEALLOCATE(convm)    ! mass flux convergence
    DEALLOCATE(w)       ! vertical velocity      
    DEALLOCATE(qv)       ! potential velocity  
    DEALLOCATE(berni)   ! bernouilli term
!!$OMP END MASTER
!!$OMP BARRIER                                                      
  END SUBROUTINE compute_caldyn
  
  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
    USE icosa
    USE vorticity_mod
    USE theta2theta_rhodz_mod
    USE pression_mod
    USE omega_mod
    USE write_field
    USE vertical_interp_mod
    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
    
    REAL(rstd) :: out_pression_lev
    CHARACTER(LEN=255) :: str_pression
    CHARACTER(LEN=255) :: physics_type
    
    out_pression_level=0
    CALL getin("out_pression_level",out_pression_level) 
    WRITE(str_pression,*) INT(out_pression_level/100)
    str_pression=ADJUSTL(str_pression)
    
    CALL writefield("ps",f_ps)
    CALL writefield("dps",f_dps)
    CALL writefield("phis",f_phis)
    CALL vorticity(f_u,f_buf_v)
    CALL writefield("vort",f_buf_v)

    CALL w_omega(f_ps, f_u, f_buf_i)
    CALL writefield('omega', f_buf_i)
    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
    ENDIF
    
    ! Temperature
    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
    
    CALL getin('physics',physics_type)
    IF (TRIM(physics_type)=='dcmip') THEN
      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
      CALL writefield("T",f_buf1_i)
      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
        CALL writefield("T"//TRIM(str_pression),f_buf_s)
      ENDIF
    ELSE
      CALL writefield("T",f_buf_i)
      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
        CALL writefield("T"//TRIM(str_pression),f_buf_s)
      ENDIF
    ENDIF
   
    ! velocity components
    CALL un2ulonlat(f_u, f_buf_i3, f_buf1_i, f_buf2_i)
    CALL writefield("ulon",f_buf1_i)
    CALL writefield("ulat",f_buf2_i)

    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
    ENDIF
    
    ! geopotential
    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
    CALL writefield("p",f_buf_p)
    CALL writefield("phi",f_buf_i)
    CALL writefield("theta",f_buf1_i) ! potential temperature
    CALL writefield("pk",f_buf2_i) ! Exner pressure
  
  
  END SUBROUTINE write_output_fields
  
  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
    USE field_mod
    USE pression_mod
    USE exner_mod
    USE geopotential_mod
    USE theta2theta_rhodz_mod
    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
         phi(:,:), phis(:), ps(:), pks(:)
    INTEGER :: ind

    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       ps = f_ps(ind)
       p  = f_p(ind)
       CALL compute_pression(ps,p,0)
       pk = f_pk(ind)
       pks = f_pks(ind)
       CALL compute_exner(ps,p,pks,pk,0)
       theta_rhodz = f_theta_rhodz(ind)
       theta = f_theta(ind)
       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
       phis = f_phis(ind)
       phi = f_phi(ind)
       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
    END DO

  END SUBROUTINE thetarhodz2geopot
  
  SUBROUTINE un2ulonlat(f_u, f_u3d, f_ulon, f_ulat)
    USE field_mod
    USE wind_mod
    TYPE(t_field), POINTER :: f_u(:), & ! IN  : normal velocity components on edges
         f_u3d(:), f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
    REAL(rstd),POINTER :: u(:,:), u3d(:,:,:), ulon(:,:), ulat(:,:)
    INTEGER :: ind
    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       u=f_u(ind)
       u3d=f_u3d(ind)
       CALL compute_wind_centered(u,u3d)
       ulon=f_ulon(ind)
       ulat=f_ulat(ind)
       CALL compute_wind_centered_lonlat_compound(u3d, ulon, ulat)
    END DO
  END SUBROUTINE un2ulonlat

  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_TV(:)
    TYPE(t_field), POINTER :: f_q(:)
    TYPE(t_field), POINTER :: f_T(:)
    
    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
    INTEGER :: ind
    
    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       Tv=f_Tv(ind)
       q=f_q(ind)
       T=f_T(ind)
       T=Tv/(1+0.608*q(:,:,1))
    END DO
    
  END SUBROUTINE Tv2T
  
END MODULE caldyn_gcm_mod