MODULE vorticity_mod

CONTAINS

  SUBROUTINE vorticity(f_ue,f_vort)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_ue(:)
    TYPE(t_field), POINTER :: f_vort(:)
  
    REAL(rstd), POINTER :: ue(:,:)
    REAL(rstd), POINTER :: vort(:,:)
    INTEGER :: ind

    CALL transfert_request(f_ue,req_e1_vect)
   
    DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      ue=f_ue(ind)
      vort=f_vort(ind)
      CALL compute_vorticity(ue, vort)
    ENDDO
  
  END SUBROUTINE vorticity
  
  SUBROUTINE compute_vorticity(ue,vort)
  USE icosa
  USE disvert_mod
  USE omp_para
  IMPLICIT NONE
    REAL(rstd),INTENT(IN)  :: ue(3*iim*jjm,llm)
    REAL(rstd),INTENT(OUT) :: vort(2*iim*jjm,llm)
    INTEGER :: i,j,ij,l

    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
        
           vort(ij+z_up,l) = 1./Av(ij+z_up)*(  ne(ij,rup)        * ue(ij+u_rup,l)        * de(ij+u_rup)         &
                                + ne(ij+t_rup,left) * ue(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)               &
                                - ne(ij,lup)        * ue(ij+u_lup,l)        * de(ij+u_lup) )                               

           vort(ij+z_down,l) = 1./Av(ij+z_down)*(  ne(ij,ldown)         * ue(ij+u_ldown,l)         * de(ij+u_ldown)          &
                                   + ne(ij+t_ldown,right) * ue(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right)                &
                                   - ne(ij,rdown)         * ue(ij+u_rdown,l)          * de(ij+u_rdown) )

        ENDDO
      ENDDO
    ENDDO
    
  END SUBROUTINE compute_vorticity

END MODULE vorticity_mod