MODULE icolmdz_param_gravity_wave
  
   INTEGER, PARAMETER :: gw_legacy=1
   INTEGER, PARAMETER :: gw_sso=2
   INTEGER,SAVE       :: gw_method
   !$OMP THREADPRIVATE(gw_method)

CONTAINS
  
  SUBROUTINE init_param_gravity_wave
  USE getin_mod
  USE xios_mod
  IMPLICIT NONE
    CHARACTER(LEN=255) :: param_gw_method
    LOGICAL :: create_etat0_limit

    create_etat0_limit = .FALSE.
    CALL getin('create_etat0_limit', create_etat0_limit)
    IF (.NOT. create_etat0_limit) RETURN

    param_gw_method='legacy'
    CALL getin('param_gw_method', param_gw_method)
    SELECT CASE (TRIM(param_gw_method))
      CASE ('legacy')
        gw_method = gw_legacy
      CASE ('sso')
        gw_method = gw_sso
      CASE default
         PRINT *, 'Bad selector for param_gw_method : <', TRIM(param_gw_method), &
            ' > options are <legacy>, <sso>'
         STOP
    END SELECT    

    SELECT CASE (gw_method)
      CASE (gw_legacy)
      !$OMP MASTER
        CALL xios_set_file_attr("relief_gw",enabled=.TRUE.)
        CALL xios_set_fieldgroup_attr("gw_read_access",read_access=.TRUE.)
      !$OMP END MASTER
      CASE (gw_sso)
      !$OMP MASTER
        CALL xios_set_file_attr("orography",enabled=.TRUE.)
        CALL xios_set_fieldgroup_attr("gwsso_read_access",read_access=.TRUE.)
      !$OMP END MASTER
    END SELECT  

  END SUBROUTINE init_param_gravity_wave

  SUBROUTINE param_gravity_wave
  USE getin_mod
  IMPLICIT NONE
   LOGICAL :: create_etat0_limit

    create_etat0_limit = .FALSE.
    CALL getin('create_etat0_limit', create_etat0_limit)
    IF (.NOT. create_etat0_limit) RETURN
      
    SELECT CASE (gw_method)
      CASE (gw_legacy)
        CALL param_gravity_wave_legacy
      CASE (gw_sso)
        CALL param_gravity_wave_sso
    END SELECT

  END SUBROUTINE param_gravity_wave

  SUBROUTINE param_gravity_wave_sso
    ! from dynamico
  USE icosa
  USE xios_mod
  USE mpipara
  USE earth_const
  USE transfert_mod
  ! from lmdz physic
  USE create_etat0_unstruct_mod, ONLY : init_param_gw_phy => init_param_gw
  USE dimphy, ONLY : klon
  !from icosa_lmdz
  USE distrib_icosa_lmdz_mod
  IMPLICIT NONE
    TYPE(t_field), POINTER, SAVE :: f_zmea(:), f_zpic(:), f_zval(:), f_zstd(:), f_zsig(:), f_zgam(:), f_zthe(:) 
     REAL, ALLOCATABLE :: zmea_phy(:), zpic_phy(:), zval_phy(:), zstd_phy(:), zsig_phy(:), zgam_phy(:), zthe_phy(:) 

    CALL allocate_field(f_zmea, field_t, type_real)
    CALL allocate_field(f_zpic, field_t, type_real)
    CALL allocate_field(f_zval, field_t, type_real)
    CALL allocate_field(f_zstd, field_t, type_real)
    CALL allocate_field(f_zsig, field_t, type_real)
    CALL allocate_field(f_zgam, field_t, type_real)
    CALL allocate_field(f_zthe, field_t, type_real)

    CALL xios_read_field("zmea_sso",f_zmea)
    CALL transfert_request(f_zmea, req_i0)
    CALL transfert_request(f_zmea, req_i1)

    CALL xios_read_field("zstd_sso",f_zstd)
    CALL transfert_request(f_zstd, req_i0)
    CALL transfert_request(f_zstd, req_i1)

    CALL xios_read_field("zsig_sso",f_zsig)
    CALL transfert_request(f_zsig, req_i0)
    CALL transfert_request(f_zsig, req_i1)

    CALL xios_read_field("zgam_sso",f_zgam)
    CALL transfert_request(f_zgam, req_i0)
    CALL transfert_request(f_zgam, req_i1)

    CALL xios_read_field("zthe_sso",f_zthe)
    CALL transfert_request(f_zthe, req_i0)
    CALL transfert_request(f_zthe, req_i1)

    ALLOCATE(zmea_phy(klon), zpic_phy(klon), zval_phy(klon), zstd_phy(klon), zsig_phy(klon), zgam_phy(klon), zthe_phy(klon))
    CALL transfer_icosa_to_lmdz(f_zmea  , zmea_phy)
    CALL transfer_icosa_to_lmdz(f_zstd  , zstd_phy)
    CALL transfer_icosa_to_lmdz(f_zsig  , zsig_phy)
    CALL transfer_icosa_to_lmdz(f_zgam  , zgam_phy)
    CALL transfer_icosa_to_lmdz(f_zthe  , zthe_phy)
    zpic_phy(:) = zmea_phy(:) + zstd_phy(:)
    zval_phy(:) = MAX(zmea_phy(:) - zstd_phy(:), 0.) 
    CALL init_param_gw_phy(zmea_phy, zpic_phy, zval_phy, zstd_phy, zsig_phy, zgam_phy, zthe_phy)

  END SUBROUTINE param_gravity_wave_sso

  SUBROUTINE param_gravity_wave_legacy
  ! from dynamico
  USE icosa
  USE xios_mod
  USE mpipara
  USE earth_const
  USE transfert_mod
  ! from lmdz physic
  USE create_etat0_unstruct_mod, ONLY : init_param_gw_phy => init_param_gw
  USE dimphy, ONLY : klon
  !from icosa_lmdz
  USE distrib_icosa_lmdz_mod
  IMPLICIT NONE
    INTEGER :: ibegin, jbegin, ni, nj, ni_glo,nj_glo
    REAL, ALLOCATABLE :: z(:,:)
    REAL, ALLOCATABLE :: lon1d(:), lon2d_in(:,:), lon2d_out(:,:)
    REAL, ALLOCATABLE :: lat1d(:), lat2d_in(:,:), lat2d_out(:,:)
    REAL,ALLOCATABLE :: delta_x(:), delta_y(:)
    REAL, ALLOCATABLE :: mask_in(:,:)
    REAL, ALLOCATABLE :: zmea_in(:,:), zpic_in(:,:), zval_in(:,:), ztz_in(:,:), zytzy_in(:,:), zxtzx_in(:,:), zxtzy_in(:,:)

    TYPE(t_field), POINTER, SAVE :: f_mask(:), f_zmea(:), f_zpic(:), f_zval(:), f_ztz(:), f_zytzy(:), f_zxtzx(:), f_zxtzy(:)
    TYPE(t_field), POINTER, SAVE :: f_zstd(:), f_zsig(:), f_zgam(:), f_zthe(:) 
    REAL(rstd),POINTER :: mask(:), zmea(:), zpic(:), zval(:), ztz(:), zytzy(:), zxtzx(:), zxtzy(:), zstd(:), zsig(:), zgam(:), zthe(:) 
    REAL, ALLOCATABLE :: zmea_phy(:), zpic_phy(:), zval_phy(:), zstd_phy(:), zsig_phy(:), zgam_phy(:), zthe_phy(:) 

    INTEGER :: i, j, ind

    CALL allocate_field(f_mask, field_t, type_real)
    CALL allocate_field(f_zmea, field_t, type_real)
    CALL allocate_field(f_zpic, field_t, type_real)
    CALL allocate_field(f_zval, field_t, type_real)
    CALL allocate_field(f_ztz, field_t, type_real)
    CALL allocate_field(f_zytzy, field_t, type_real)
    CALL allocate_field(f_zxtzx, field_t, type_real)
    CALL allocate_field(f_zxtzy, field_t, type_real)
    CALL allocate_field(f_zstd, field_t, type_real)
    CALL allocate_field(f_zsig, field_t, type_real)
    CALL allocate_field(f_zgam, field_t, type_real)
    CALL allocate_field(f_zthe, field_t, type_real)

!$OMP MASTER
    CALL xios_get_domain_attr("domain_relief_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj,  nj_glo=nj_glo)
    PRINT*, ni,nj,ni_glo,nj_glo
    
    ALLOCATE( lon1d(ni), lat1d(nj))
    CALL xios_get_domain_attr("domain_relief_gw",lonvalue_1d=lon1d, latvalue_1d=lat1d)
       
    ALLOCATE(z(0:ni+1,0:nj+1))
    CALL xios_recv_field("relief_exp", z)

    ALLOCATE(lon2d_in(ni,nj), lat2d_in(ni,nj) )
    ALLOCATE(lon2d_out(0:ni+1,0:nj+1), lat2d_out(0:ni+1,0:nj+1) )
    
    DO i=1,ni
     lon2d_in(i,:) = lon1d(i)
    END DO

    DO j=1,nj
     lat2d_in(:,j) = lat1d(j)
    END DO

    CALL xios_send_field("lon2d_in", lon2d_in)
    CALL xios_recv_field("lon2d_out", lon2d_out)
    CALL xios_send_field("lat2d_in", lat2d_in)
    CALL xios_recv_field("lat2d_out", lat2d_out)
   
    ALLOCATE(delta_x(nj))
    ALLOCATE(delta_y(nj))
    
    delta_y(:) = Pi*radius / nj_glo  !! difference with  lmdz reference : delta_y=2*Pi*radius/nj_glo ??
    
    DO j=1,nj
      delta_x(j) = 2*Pi*radius / ni_glo * cos( lat2d_out(1,j) * pi / 180.)
    ENDDO
    
    !north pole
    IF (jbegin==0) THEN
      z(:,0) = z(:,1)
      delta_y(1) = delta_y(1)*0.5
    ENDIF
    
    !south pole
    IF (jbegin+nj==nj_glo) THEN
      z(:,nj+1) = z(:,nj)
      delta_y(nj) = delta_y(nj)*0.5
    ENDIF
    
    ALLOCATE(mask_in(ni,nj))
    ALLOCATE(zmea_in(ni,nj), zpic_in(ni,nj), zval_in(ni,nj), ztz_in(ni,nj), zytzy_in(ni,nj), zxtzx_in(ni,nj), zxtzy_in(ni,nj) )
    
    mask_in(:,:)=0
    DO j=1,nj
      DO i=1,ni
        IF ( z(i,j)>1. ) mask_in(i,j) = 1.
        zmea_in(i,j)  = z(i,j)
        zpic_in(i,j)  = z(i,j)
        zval_in(i,j)  = z(i,j)
        ztz_in(i,j)   = z(i,j)*z(i,j)
        zytzy_in(i,j) = (z(i,j+1)-z(i,j-1))**2/(2*delta_y(j))**2
        zxtzx_in(i,j) = (z(i+1,j)-z(i-1,j))**2/(2*delta_x(j))**2
        zxtzy_in(i,j) = (z(i,j+1)-z(i,j-1)) / (2*delta_y(j)) *(z(i+1,j)-z(i-1,j)) / (2*delta_x(j))
      ENDDO
    ENDDO

    CALL xios_send_field("mask_in",mask_in)
    CALL xios_send_field("zmea_in",zmea_in)
    CALL xios_send_field("zpic_in",zpic_in)
    CALL xios_send_field("zval_in",zval_in)
    CALL xios_send_field("ztz_in",ztz_in)
    CALL xios_send_field("zytzy_in",zytzy_in)
    CALL xios_send_field("zxtzx_in",zxtzx_in)
    CALL xios_send_field("zxtzy_in",zxtzy_in)
    
    !! only for test, to compare projection on ico grid Vs lmdz_reg grid
    !! CALL compute_regular_param

    !$OMP END MASTER
    
    CALL xios_read_field("mask_out", f_mask)
    CALL transfert_request(f_mask, req_i0)
    CALL transfert_request(f_mask, req_i1)

    CALL xios_read_field("zmea_out",f_zmea)
    CALL transfert_request(f_zmea, req_i0)
    CALL transfert_request(f_zmea, req_i1)

    CALL xios_read_field("zpic_out",f_zpic)
    CALL transfert_request(f_zpic, req_i0)
    CALL transfert_request(f_zpic, req_i1)

    CALL xios_read_field("zval_out",f_zval)
    CALL transfert_request(f_zval, req_i0)
    CALL transfert_request(f_zval, req_i1)

    CALL xios_read_field("ztz_out", f_ztz)
    CALL transfert_request(f_ztz, req_i0)
    CALL transfert_request(f_ztz, req_i1)

    CALL xios_read_field("zytzy_out", f_zytzy)
    CALL transfert_request(f_zytzy, req_i0)
    CALL transfert_request(f_zytzy, req_i1)

    CALL xios_read_field("zxtzx_out", f_zxtzx)
    CALL transfert_request(f_zxtzx, req_i0)
    CALL transfert_request(f_zxtzx, req_i1)

    CALL xios_read_field("zxtzy_out", f_zxtzy)
    CALL transfert_request(f_zxtzy, req_i0)
    CALL transfert_request(f_zxtzy, req_i1)

   !$OMP BARRIER
   !$OMP MASTER

   !! compute zstd

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      zmea =f_zmea(ind)
      zstd = f_zstd(ind)
      ztz  = f_ztz(ind)
      CALL compute_zstd(zmea, ztz, zstd)
    ENDDO

   !! smmothing ?

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      
      zmea=f_zmea(ind)
      zpic=f_zpic(ind)
      zval=f_zval(ind)
      zstd=f_zstd(ind)
      zytzy=f_zytzy(ind)
      zxtzx=f_zxtzx(ind)
      zxtzy=f_zxtzy(ind)
 
      CALL compute_smoothing(zmea) 
      CALL compute_smoothing(zpic) 
      CALL compute_smoothing(zval) 
      CALL compute_smoothing(zstd) 
      CALL compute_smoothing(zytzy) 
      CALL compute_smoothing(zxtzx) 
      CALL compute_smoothing(zxtzy) 
 
    ENDDO

    !! SLOPE, ANISOTROPY AND THETA ANGLE
    
    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)

      zytzy=f_zytzy(ind)
      zxtzx=f_zxtzx(ind)
      zxtzy=f_zxtzy(ind)
      zsig=f_zsig(ind)
      zgam=f_zgam(ind)
      zthe=f_zthe(ind) 
      CALL compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe) 
   
    ENDDO

    !! apply mask
    
    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)

      mask=f_mask(ind)
      zmea=f_zmea(ind)
      zpic=f_zpic(ind)
      zval=f_zval(ind)
      zstd=f_zstd(ind)
      zsig=f_zsig(ind)
      zgam=f_zgam(ind)
      zthe=f_zthe(ind) 

      CALL compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe) 

    ENDDO

!$OMP END MASTER
!$OMP BARRIER
  
  !! for debugging
  !!  CALL xios_write_field("zmask", f_mask)
  !!  CALL xios_write_field("zmea", f_zmea)
  !!  CALL xios_write_field("zpic", f_zpic)
  !!  CALL xios_write_field("zval", f_zval)
  !!  CALL xios_write_field("zsig", f_zsig)
  !!  CALL xios_write_field("zgam", f_zgam)
  !!  CALL xios_write_field("zthe", f_zthe)
    
    ALLOCATE(zmea_phy(klon), zpic_phy(klon), zval_phy(klon), zstd_phy(klon), zsig_phy(klon), zgam_phy(klon), zthe_phy(klon))
    CALL transfer_icosa_to_lmdz(f_zmea  , zmea_phy)
    CALL transfer_icosa_to_lmdz(f_zpic  , zpic_phy)
    CALL transfer_icosa_to_lmdz(f_zval  , zval_phy)
    CALL transfer_icosa_to_lmdz(f_zstd  , zstd_phy)
    CALL transfer_icosa_to_lmdz(f_zsig  , zsig_phy)
    CALL transfer_icosa_to_lmdz(f_zgam  , zgam_phy)
    CALL transfer_icosa_to_lmdz(f_zthe  , zthe_phy)
    CALL init_param_gw_phy(zmea_phy, zpic_phy, zval_phy, zstd_phy, zsig_phy, zgam_phy, zthe_phy)

  END SUBROUTINE param_gravity_wave_legacy

  SUBROUTINE compute_smoothing(X)
  USE icosa
  IMPLICIT NONE
    REAL(rstd),INTENT(INOUT)  :: X(iim*jjm)
    REAL(rstd)                :: tmp(iim*jjm)
    INTEGER :: i,j,ij

    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        tmp(ij) = (X(ij) + (X(ij+t_right) + X(ij+t_rup) + X(ij+t_lup) + X(ij+t_left) + X(ij+t_ldown) + X(ij+t_rdown))/3.)/3. 
      ENDDO
    ENDDO
  
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        X(ij) = tmp(ij) 
      ENDDO
    ENDDO
  
  END SUBROUTINE

  SUBROUTINE compute_zstd(zmea, ztz, zstd)
  USE icosa
  IMPLICIT NONE
    REAL(rstd),INTENT(IN)   :: zmea(iim*jjm)
    REAL(rstd),INTENT(IN)   :: ztz(iim*jjm)
    REAL(rstd),INTENT(OUT) :: zstd(iim*jjm)

    INTEGER :: i,j,ij

    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        zstd(ij) = ztz(ij) - zmea(ij)*zmea(ij)
        IF (zstd(ij)<0) zstd(ij)=0
        zstd(ij)=SQRT(zstd(ij))
      ENDDO
    ENDDO
  
  END SUBROUTINE compute_zstd


  SUBROUTINE compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe)
  USE icosa
  IMPLICIT NONE
    REAL(rstd),INTENT(IN)   :: zxtzx(iim*jjm)
    REAL(rstd),INTENT(IN)   :: zytzy(iim*jjm)
    REAL(rstd),INTENT(IN)   :: zxtzy(iim*jjm)
    REAL(rstd),INTENT(OUT)  :: zsig(iim*jjm)
    REAL(rstd),INTENT(OUT)  :: zgam(iim*jjm)
    REAL(rstd),INTENT(OUT)  :: zthe(iim*jjm)
    
    REAL(rstd) :: xk,xl,xm,xp,xq,xw

    INTEGER :: i,j,ij

    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
      
        ij=(j-1)*iim+i
        xk=(zxtzx(ij)+zytzy(ij))/2.
        xl=(zxtzx(ij)-zytzy(ij))/2.
        xm=zxtzy(ij)
        xp=xk-SQRT(xl**2+xm**2)
        xq=xk+SQRT(xl**2+xm**2)
        xw=1.e-8
        IF(xp<=xw) xp=0.
        IF(xq<=xw) xq=xw
        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
        
        !--- SLOPE, ANISOTROPY AND THETA ANGLE
      
        zsig(ij)=SQRT(xq)
        zgam(ij)=xp/xq
        zthe(ij)=90.*ATAN2(xm,xl)/Pi

      ENDDO
    ENDDO

  END SUBROUTINE compute_sigma_gamma_theta

  SUBROUTINE compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe) 
  USE icosa
  IMPLICIT NONE
    REAL(rstd),INTENT(IN)   :: mask(iim*jjm)
    REAL(rstd),INTENT(OUT)   :: zmea(iim*jjm)
    REAL(rstd),INTENT(OUT)   :: zpic(iim*jjm)
    REAL(rstd),INTENT(OUT)   :: zval(iim*jjm)
    REAL(rstd),INTENT(OUT)   :: zstd(iim*jjm)
    REAL(rstd),INTENT(OUT)   :: zsig(iim*jjm)
    REAL(rstd),INTENT(OUT)   :: zgam(iim*jjm)
    REAL(rstd),INTENT(OUT)   :: zthe(iim*jjm)
    

    INTEGER :: i,j,ij

    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        IF (mask(ij)<0.1) THEN
          zmea(ij)=0.
          zpic(ij)=0.
          zval(ij)=0.
          zstd(ij)=0.
          zsig(ij)=0.
          zgam(ij)=0.
          zthe(ij)=0.
        ENDIF
      ENDDO
    ENDDO
  
  END SUBROUTINE compute_masking

  SUBROUTINE compute_regular_param
  USE xios_mod
  USE icosa
  IMPLICIT NONE
    REAL, ALLOCATABLE :: mask(:,:) , zmea(:,:), zpic(:,:), zval(:,:), ztz(:,:), zstd(:,:)
    REAL, ALLOCATABLE :: zytzy(:,:) , zxtzx(:,:), zxtzy(:,:), zsig(:,:) , zgam(:,:), zthe(:,:)
    REAL :: xk, xl, xm, xp, xq, xw
    INTEGER :: ni, nj, ibegin, jbegin, ni_glo, nj_glo
    INTEGER :: i,j

    CALL xios_get_domain_attr("regular_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj,  nj_glo=nj_glo)

    ALLOCATE(mask(0:ni+1,0:nj+1), zmea(0:ni+1,0:nj+1), zpic(0:ni+1,0:nj+1), zval(0:ni+1,0:nj+1), ztz(0:ni+1,0:nj+1))
    ALLOCATE(zstd(0:ni+1,0:nj+1), zytzy(0:ni+1,0:nj+1), zxtzx(0:ni+1,0:nj+1), zxtzy(0:ni+1,0:nj+1))
    ALLOCATE(zsig(0:ni+1,0:nj+1), zgam(0:ni+1,0:nj+1), zthe(0:ni+1,0:nj+1))

    CALL xios_recv_field("mask_reg_exp",mask)
    CALL xios_recv_field("zmea_reg_exp",zmea)
    CALL xios_recv_field("zpic_reg_exp",zpic)
    CALL xios_recv_field("zval_reg_exp",zval)
    CALL xios_recv_field("ztz_reg_exp",ztz)
    CALL xios_recv_field("zytzy_reg_exp",zytzy)
    CALL xios_recv_field("zxtzx_reg_exp",zxtzx)
    CALL xios_recv_field("zxtzy_reg_exp",zxtzy)
    
    DO j=0,nj+1
      DO i=0, ni+1
        zstd(i,j) = ztz(i,j)-zmea(i,j)*zmea(i,j)
      ENDDO
    ENDDO

    ! smoothing
    CALL smoothing_regular(ni,nj, zmea)
    CALL smoothing_regular(ni,nj, zpic)
    CALL smoothing_regular(ni,nj, zval)
    CALL smoothing_regular(ni,nj, zstd)
    CALL smoothing_regular(ni,nj, zytzy)
    CALL smoothing_regular(ni,nj, zxtzx)
    CALL smoothing_regular(ni,nj, zxtzy)

    DO j=1,nj
      DO i=1, ni
        xk=(zxtzx(i,j)+zytzy(i,j))/2.
        xl=(zxtzx(i,j)-zytzy(i,j))/2.
        xm=zxtzy(i,j)
        xp=xk-SQRT(xl**2+xm**2)
        xq=xk+SQRT(xl**2+xm**2)
        xw=1.e-8
        IF(xp<=xw) xp=0.
        IF(xq<=xw) xq=xw
        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
        
        !--- SLOPE, ANISOTROPY AND THETA ANGLE
      
        zsig(i,j)=SQRT(xq)
        zgam(i,j)=xp/xq
        zthe(i,j)=90.*ATAN2(xm,xl)/Pi
      ENDDO
    ENDDO

    DO j=1,nj
      DO i=1,ni
        IF (mask(i,j)<0.1) THEN
          zmea(i,j)=0.
          zpic(i,j)=0.
          zval(i,j)=0.
          zstd(i,j)=0.
          zsig(i,j)=0.
          zgam(i,j)=0.
          zthe(i,j)=0.
        ENDIF
      ENDDO
    ENDDO

    CALL xios_send_field("zmask_regout", mask(1:ni,1:nj))
    CALL xios_send_field("zmea_regout", zmea(1:ni,1:nj))
    CALL xios_send_field("zpic_regout", zpic(1:ni,1:nj))
    CALL xios_send_field("zval_regout", zval(1:ni,1:nj))
    CALL xios_send_field("zsig_regout", zsig(1:ni,1:nj))
    CALL xios_send_field("zgam_regout", zgam(1:ni,1:nj))
    CALL xios_send_field("zthe_regout", zthe(1:ni,1:nj))

  END SUBROUTINE compute_regular_param

  SUBROUTINE smoothing_regular(ni,nj,var)
  IMPLICIT NONE
    INTEGER, INTENT(IN) :: ni, nj
    REAL, INTENT(INOUT) :: var(0:ni+1,0:nj+1)
    REAL :: tmp(0:ni+1,0:nj+1)
    INTEGER :: i,j

    DO j=1,nj
      DO i=1,ni
        tmp(i,j) = (var(i,j) + 0.5*(var(i+1,j) + var(i-1,j) + var(i,j+1) +var(i,j-1))              &
                             + 0.25*(var(i+1,j+1) + var(i-1,j+1) + var(i+1,j-1) + var(i-1,j-1)))/4.
      ENDDO
    ENDDO
    
    DO j=1,nj
      DO i=1,ni
        var(i,j)=tmp(i,j)
      ENDDO
    ENDDO

  END SUBROUTINE smoothing_regular


END MODULE icolmdz_param_gravity_wave