Commit 3de69958 authored by dubos's avatar dubos
Browse files

devel : diagnose divergence and vorticity

git-svn-id: svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/dynamico/svn/codes/icosagcm/devel@1052 5a24c8df-7893-42df-8d02-8693bdefd653
parent 5778cc8f
# module c++/gnu is needed to link with XIOS
# module feature/openmpi/mpi_compiler/intel is made necessary by c++/gnu
module purge
module load feature/openmpi/mpi_compiler/intel
module load c++/gnu/7.3.0
module load intel/17.0.6.256
module load mpi/openmpi/2.0.4
module unload netcdf-c netcdf-fortran hdf5 flavor perl hdf5 boost blitz mpi gnu
#module purge
module load flavor/hdf5/parallel
module load netcdf-fortran/4.4.4
module load mkl/17.0.6.256
module load hdf5/1.8.20
# module gnu is needed to link with XIOS
# module feature/openmpi/mpi_compiler/intel is made necessary by gnu
module load gnu
module load feature/bridge/heterogenous_mpmd
module load parmetis/4.0.3
......@@ -54,6 +54,8 @@ CONTAINS
USE compute_rhodz_mod
USE compute_pression_mod
USE compute_temperature_mod
USE compute_vorticity_mod
USE compute_divergence_mod
USE vertical_interp_mod
USE compute_caldyn_mod
......@@ -67,6 +69,8 @@ CONTAINS
compute_pression => compute_pression_hex
compute_pression_mid => compute_pression_mid_hex
compute_temperature => compute_temperature_hex
compute_vorticity => compute_vorticity_hex
compute_divergence => compute_divergence_hex
compute_hydrostatic_pressure => compute_hydrostatic_pressure_hex
compute_vertical_interp => compute_vertical_interp_hex
compute_pvort_only => compute_pvort_only_hex
......@@ -82,6 +86,8 @@ CONTAINS
USE compute_rhodz_mod
USE compute_pression_mod
USE compute_temperature_mod
USE compute_vorticity_mod
USE compute_divergence_mod
USE vertical_interp_mod
USE compute_caldyn_mod
......@@ -95,6 +101,8 @@ CONTAINS
compute_pression => compute_pression_unst
compute_pression_mid => compute_pression_mid_unst
compute_temperature => compute_temperature_unst
compute_vorticity => compute_vorticity_unst
compute_divergence => compute_divergence_unst
compute_hydrostatic_pressure => compute_hydrostatic_pressure_unst
compute_vertical_interp => compute_vertical_interp_unst
compute_pvort_only => compute_pvort_only_unst
......
......@@ -49,7 +49,7 @@ CONTAINS
END SUBROUTINE init_grid_param
{%- set diags=('rhodz','pression','pression_mid','temperature','hydrostatic_pressure','vertical_interp',
{%- set diags=('rhodz','pression','pression_mid','temperature','vorticity','divergence', 'hydrostatic_pressure','vertical_interp',
'pvort_only','theta','geopot','caldyn_fast','caldyn_slow_hydro','caldyn_coriolis') %}
{% for grid in 'hex','unst' %}
SUBROUTINE select_compute_{{grid}}
......@@ -57,6 +57,8 @@ CONTAINS
USE compute_rhodz_mod
USE compute_pression_mod
USE compute_temperature_mod
USE compute_vorticity_mod
USE compute_divergence_mod
USE vertical_interp_mod
USE compute_caldyn_mod
......
......@@ -51,7 +51,7 @@ CONTAINS
!--------------------------------- Basic check --------------------------------
SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it)
USE vorticity_mod
USE compute_vorticity_mod
USE caldyn_gcm_mod
USE exner_mod
USE omp_para, ONLY : is_master
......@@ -152,10 +152,8 @@ CONTAINS
!------------------------ Detailed check (Lauritzen et al., 2014) -----------------------------
SUBROUTINE check_conserve_detailed(it,tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis)
USE vorticity_mod
USE caldyn_gcm_mod
USE exner_mod
! USE mpipara, ONLY : is_mpi_root, comm_icosa
USE omp_para, ONLY : is_master
INTEGER, INTENT(IN) :: tag
TYPE(t_field),POINTER :: f_ps(:)
......@@ -267,7 +265,6 @@ CONTAINS
SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot)
USE vorticity_mod
TYPE(t_field), POINTER :: f_ue(:)
TYPE(t_field), POINTER :: f_theta_rhodz(:)
TYPE(t_field), POINTER :: f_phis(:)
......
......@@ -26,6 +26,18 @@ MODULE compute_diagnostics_mod
REAL(rstd),INTENT(INOUT) :: temp(:,:)
END SUBROUTINE comp_temperature
SUBROUTINE comp_vorticity(ue, vort)
IMPORT
REAL(rstd),INTENT(IN) :: ue(:,:)
REAL(rstd),INTENT(OUT) :: vort(:,:)
END SUBROUTINE comp_vorticity
SUBROUTINE comp_divergence(ue, div)
IMPORT
REAL(rstd),INTENT(IN) :: ue(:,:)
REAL(rstd),INTENT(OUT) :: div(:,:)
END SUBROUTINE comp_divergence
SUBROUTINE comp_hydro_press(rhodz, theta_rhodz, ps, p)
IMPORT
REAL(rstd),INTENT(IN) :: rhodz(:,:)
......@@ -47,10 +59,12 @@ MODULE compute_diagnostics_mod
PROCEDURE(comp_rhodz), POINTER, SAVE :: compute_rhodz => NULL()
PROCEDURE(comp_pression), POINTER, SAVE :: compute_pression => NULL(), compute_pression_mid => NULL()
PROCEDURE(comp_temperature), POINTER, SAVE :: compute_temperature => NULL()
PROCEDURE(comp_vorticity), POINTER, SAVE :: compute_vorticity => NULL()
PROCEDURE(comp_divergence), POINTER, SAVE :: compute_divergence => NULL()
PROCEDURE(comp_hydro_press), POINTER, SAVE :: compute_hydrostatic_pressure => NULL()
PROCEDURE(comp_vert_interp), POINTER, SAVE :: compute_vertical_interp => NULL()
PUBLIC :: compute_rhodz, compute_pression, compute_pression_mid, compute_temperature, compute_hydrostatic_pressure, &
compute_vertical_interp
PUBLIC :: compute_rhodz, compute_pression, compute_pression_mid, compute_temperature, &
compute_vorticity, compute_divergence, compute_hydrostatic_pressure, compute_vertical_interp
END MODULE compute_diagnostics_mod
MODULE compute_divergence_mod
USE grid_param
USE prec, ONLY : rstd
IMPLICIT NONE
PRIVATE
#include "../unstructured/unstructured.h90"
PUBLIC :: divergence, compute_divergence_hex, compute_divergence_unst
CONTAINS
#ifdef BEGIN_DYSL
KERNEL(compute_divergence)
FORALL_CELLS_EXT()
ON_PRIMAL
div_ij = 0.d0
FORALL_EDGES
div_ij = div_ij + SIGN*ue(EDGE)*LE
END_BLOCK
div(CELL) = div_ij / AI
END_BLOCK
END_BLOCK
END_BLOCK
#endif END_DYSL
SUBROUTINE divergence(f_ue,f_div)
USE icosa
USE compute_diagnostics_mod, ONLY : compute_divergence
TYPE(t_field), POINTER :: f_ue(:)
TYPE(t_field), POINTER :: f_div(:)
REAL(rstd), POINTER :: ue(:,:)
REAL(rstd), POINTER :: div(:,:)
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)
div=f_div(ind)
CALL compute_divergence(ue, div)
ENDDO
END SUBROUTINE divergence
!-------------- Wrappers for F2008 conformity -----------------
SUBROUTINE compute_divergence_unst(ue,div)
REAL(rstd),INTENT(IN) :: ue(:,:)
REAL(rstd),INTENT(OUT) :: div(:,:)
CALL compute_divergence_unst_(ue,div)
END SUBROUTINE compute_divergence_unst
SUBROUTINE compute_divergence_hex(ue,div)
REAL(rstd),INTENT(IN) :: ue(:,:)
REAL(rstd),INTENT(OUT) :: div(:,:)
CALL compute_divergence_hex_(ue,div)
END SUBROUTINE compute_divergence_hex
!--------------------------------------------------------------
SUBROUTINE compute_divergence_unst_(ue,div)
USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
USE geometry, ONLY : le, Ai
USE data_unstructured_mod, ONLY : primal_deg, primal_edge, primal_ne
FIELD_U :: ue
FIELD_MASS :: div
DECLARE_INDICES
DECLARE_EDGES
DECLARE_VERTICES
NUM :: div_ij
#include "../kernels_unst/compute_divergence.k90"
END SUBROUTINE compute_divergence_unst_
SUBROUTINE compute_divergence_hex_(ue,div)
USE icosa
USE omp_para
REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
REAL(rstd),INTENT(OUT) :: div(iim*jjm,llm)
REAL(rstd) :: div_ij
INTEGER :: ij,l
#include "../kernels_hex/compute_divergence.k90"
END SUBROUTINE compute_divergence_hex_
END MODULE compute_divergence_mod
MODULE compute_vorticity_mod
USE grid_param
USE prec, ONLY : rstd
IMPLICIT NONE
PRIVATE
#include "../unstructured/unstructured.h90"
PUBLIC :: vorticity, compute_vorticity_manual, compute_vorticity_hex, compute_vorticity_unst
CONTAINS
#ifdef BEGIN_DYSL
KERNEL(compute_vorticity)
FORALL_CELLS_EXT()
ON_DUAL
etav = 0.d0
FORALL_EDGES
etav = etav + SIGN*ue(EDGE)*DE
END_BLOCK
vort(DUAL_CELL) = etav / AV
END_BLOCK
END_BLOCK
END_BLOCK
#endif END_DYSL
SUBROUTINE vorticity(f_ue,f_vort)
USE icosa
USE compute_diagnostics_mod, ONLY : compute_vorticity
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
!-------------- Wrappers for F2008 conformity -----------------
SUBROUTINE compute_vorticity_unst(ue,vort)
REAL(rstd),INTENT(IN) :: ue(:,:)
REAL(rstd),INTENT(OUT) :: vort(:,:)
CALL compute_vorticity_unst_(ue,vort)
END SUBROUTINE compute_vorticity_unst
SUBROUTINE compute_vorticity_hex(ue,vort)
REAL(rstd),INTENT(IN) :: ue(:,:)
REAL(rstd),INTENT(OUT) :: vort(:,:)
CALL compute_vorticity_hex_(ue,vort)
END SUBROUTINE compute_vorticity_hex
!--------------------------------------------------------------
SUBROUTINE compute_vorticity_unst_(ue,vort)
USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
USE geometry, ONLY : de, Av
USE data_unstructured_mod, ONLY : dual_deg, dual_edge, dual_ne
FIELD_U :: ue
FIELD_Z :: vort
DECLARE_INDICES
DECLARE_EDGES
DECLARE_VERTICES
NUM :: etav
#include "../kernels_unst/compute_vorticity.k90"
END SUBROUTINE compute_vorticity_unst_
SUBROUTINE compute_vorticity_hex_(ue,vort)
USE icosa
USE omp_para
REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
REAL(rstd),INTENT(OUT) :: vort(2*iim*jjm,llm)
REAL(rstd) :: etav
INTEGER :: ij,l
#include "../kernels_hex/compute_vorticity.k90"
END SUBROUTINE compute_vorticity_hex_
SUBROUTINE compute_vorticity_manual(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_manual
END MODULE compute_vorticity_mod
......@@ -45,6 +45,8 @@ CONTAINS
USE compute_pression_mod, ONLY : pression_mid, hydrostatic_pressure
USE compute_temperature_mod
USE compute_velocity_mod
USE compute_vorticity_mod
USE compute_divergence_mod
USE vertical_interp_mod
USE theta2theta_rhodz_mod
USE compute_omega_mod, ONLY : w_omega
......@@ -108,6 +110,7 @@ CONTAINS
CALL output_field("q_init",f_q)
CALL output_field("temp_init",f_buf_i)
CALL output_field("vort_init",f_buf_v)
ELSE
CALL output_field("p",f_pmid)
CALL output_field("ps",f_ps)
......@@ -123,6 +126,16 @@ CONTAINS
CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,preff)
CALL output_field("SST",f_buf_s)
END IF
CALL vorticity(f_u, f_buf_v)
CALL divergence(f_u, f_buf_i)
IF(init) THEN
CALL output_field("vort_init",f_buf_v)
CALL output_field("div_init",f_buf_i)
ELSE
CALL output_field("vort",f_buf_v)
CALL output_field("div",f_buf_i)
END IF
IF(grid_type == grid_unst) RETURN
......
!--------------------------------------------------------------------------
!---------------------------- compute_divergence ----------------------------------
DO l = ll_begin, ll_end
!DIR$ SIMD
DO ij=ij_begin_ext, ij_end_ext
div_ij = 0.d0
div_ij = div_ij + ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
div_ij = div_ij + ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
div_ij = div_ij + ne_left*ue(ij+u_left,l)*le(ij+u_left)
div_ij = div_ij + ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
div_ij = div_ij + ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
div_ij = div_ij + ne_right*ue(ij+u_right,l)*le(ij+u_right)
div(ij,l) = div_ij / Ai(ij)
END DO
END DO
!---------------------------- compute_divergence ----------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!---------------------------- compute_vorticity ----------------------------------
DO l = ll_begin, ll_end
!DIR$ SIMD
DO ij=ij_begin_ext, ij_end_ext
etav = 0.d0
etav = etav + ne_rup*ue(ij+u_rup,l)*de(ij+u_rup)
etav = etav + ne_left*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)
etav = etav + (-ne_lup)*ue(ij+u_lup,l)*de(ij+u_lup)
vort(ij+z_up,l) = etav / Av(ij+z_up)
etav = 0.d0
etav = etav + (-ne_rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown)
etav = etav + ne_right*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)
etav = etav + ne_ldown*ue(ij+u_ldown,l)*de(ij+u_ldown)
vort(ij+z_down,l) = etav / Av(ij+z_down)
END DO
END DO
!---------------------------- compute_vorticity ----------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!---------------------------- compute_divergence ----------------------------------
DO l = 1, llm
!DIR$ SIMD
DO ij=ij_begin_ext, ij_end_ext
div_ij = 0.d0
div_ij = div_ij + ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
div_ij = div_ij + ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
div_ij = div_ij + ne_left*ue(ij+u_left,l)*le(ij+u_left)
div_ij = div_ij + ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
div_ij = div_ij + ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
div_ij = div_ij + ne_right*ue(ij+u_right,l)*le(ij+u_right)
div(ij,l) = div_ij / Ai(ij)
END DO
END DO
!---------------------------- compute_divergence ----------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!---------------------------- compute_vorticity ----------------------------------
DO l = 1, llm
!DIR$ SIMD
DO ij=ij_begin_ext, ij_end_ext
etav = 0.d0
etav = etav + ne_rup*ue(ij+u_rup,l)*de(ij+u_rup)
etav = etav + ne_left*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)
etav = etav + (-ne_lup)*ue(ij+u_lup,l)*de(ij+u_lup)
vort(ij+z_up,l) = etav / Av(ij+z_up)
etav = 0.d0
etav = etav + (-ne_rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown)
etav = etav + ne_right*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)
etav = etav + ne_ldown*ue(ij+u_ldown,l)*de(ij+u_ldown)
vort(ij+z_down,l) = etav / Av(ij+z_down)
END DO
END DO
!---------------------------- compute_vorticity ----------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!---------------------------- compute_divergence ----------------------------------
!$OMP DO SCHEDULE(STATIC)
DO ij = 1, primal_num
! this VLOOP iterates over primal cell edges
SELECT CASE(primal_deg(ij))
CASE(4)
edge1 = primal_edge(1,ij)
edge2 = primal_edge(2,ij)
edge3 = primal_edge(3,ij)
edge4 = primal_edge(4,ij)
le1 = le(edge1)
le2 = le(edge2)
le3 = le(edge3)
le4 = le(edge4)
sign1 = primal_ne(1,ij)
sign2 = primal_ne(2,ij)
sign3 = primal_ne(3,ij)
sign4 = primal_ne(4,ij)
!DIR$ SIMD
DO l = 1, llm
div_ij = 0.d0
div_ij = div_ij + sign1*ue(l,edge1)*le1
div_ij = div_ij + sign2*ue(l,edge2)*le2
div_ij = div_ij + sign3*ue(l,edge3)*le3
div_ij = div_ij + sign4*ue(l,edge4)*le4
div(l,ij) = div_ij / Ai(ij)
END DO
CASE(5)
edge1 = primal_edge(1,ij)
edge2 = primal_edge(2,ij)
edge3 = primal_edge(3,ij)
edge4 = primal_edge(4,ij)
edge5 = primal_edge(5,ij)
le1 = le(edge1)
le2 = le(edge2)
le3 = le(edge3)
le4 = le(edge4)
le5 = le(edge5)
sign1 = primal_ne(1,ij)
sign2 = primal_ne(2,ij)
sign3 = primal_ne(3,ij)
sign4 = primal_ne(4,ij)
sign5 = primal_ne(5,ij)
!DIR$ SIMD
DO l = 1, llm
div_ij = 0.d0
div_ij = div_ij + sign1*ue(l,edge1)*le1
div_ij = div_ij + sign2*ue(l,edge2)*le2
div_ij = div_ij + sign3*ue(l,edge3)*le3
div_ij = div_ij + sign4*ue(l,edge4)*le4
div_ij = div_ij + sign5*ue(l,edge5)*le5
div(l,ij) = div_ij / Ai(ij)
END DO
CASE(6)
edge1 = primal_edge(1,ij)
edge2 = primal_edge(2,ij)
edge3 = primal_edge(3,ij)
edge4 = primal_edge(4,ij)
edge5 = primal_edge(5,ij)
edge6 = primal_edge(6,ij)
le1 = le(edge1)
le2 = le(edge2)
le3 = le(edge3)
le4 = le(edge4)
le5 = le(edge5)
le6 = le(edge6)
sign1 = primal_ne(1,ij)
sign2 = primal_ne(2,ij)
sign3 = primal_ne(3,ij)
sign4 = primal_ne(4,ij)
sign5 = primal_ne(5,ij)
sign6 = primal_ne(6,ij)
!DIR$ SIMD
DO l = 1, llm
div_ij = 0.d0
div_ij = div_ij + sign1*ue(l,edge1)*le1
div_ij = div_ij + sign2*ue(l,edge2)*le2
div_ij = div_ij + sign3*ue(l,edge3)*le3
div_ij = div_ij + sign4*ue(l,edge4)*le4
div_ij = div_ij + sign5*ue(l,edge5)*le5
div_ij = div_ij + sign6*ue(l,edge6)*le6
div(l,ij) = div_ij / Ai(ij)
END DO
CASE DEFAULT
!DIR$ SIMD
DO l = 1, llm
div_ij = 0.d0
DO iedge = 1, primal_deg(ij)
edge = primal_edge(iedge,ij)
div_ij = div_ij + primal_ne(iedge,ij)*ue(l,edge)*le(edge)
END DO
div(l,ij) = div_ij / Ai(ij)
END DO
END SELECT
END DO
!$OMP END DO
!---------------------------- compute_divergence ----------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!---------------------------- compute_vorticity ----------------------------------
!$OMP DO SCHEDULE(STATIC)
DO ij = 1, dual_num
! this VLOOP iterates over dual cell edges
SELECT CASE(dual_deg(ij))
CASE(3)
edge1 = dual_edge(1,ij)
edge2 = dual_edge(2,ij)
edge3 = dual_edge(3,ij)
de1 = de(edge1)
de2 = de(edge2)
de3 = de(edge3)
sign1 = dual_ne(1,ij)
sign2 = dual_ne(2,ij)
sign3 = dual_ne(3,ij)
!DIR$ SIMD
DO l = 1, llm
etav = 0.d0
etav = etav + sign1*ue(l,edge1)*de1
etav = etav + sign2*ue(l,edge2)*de2
etav = etav + sign3*ue(l,edge3)*de3
vort(l,ij) = etav / Av(ij)
END DO
CASE(4)
edge1 = dual_edge(1,ij)
edge2 = dual_edge(2,ij)
edge3 = dual_edge(3,ij)
edge4 = dual_edge(4,ij)
de1 = de(edge1)
de2 = de(edge2)