source: codes/icosagcm/devel/src/diagnostics/diagflux.F90 @ 595

MODULE diagflux_mod
  USE icosa
  USE omp_para
  IMPLICIT NONE
  SAVE
  PRIVATE

  TYPE(t_field), POINTER, PUBLIC :: &
       f_masst(:), f_qmasst(:), & ! time-averaged mass, tracer mass, 
       f_massfluxt(:), f_qfluxt(:), & ! time-integrated mass flux and tracer flux
       f_qfluxt_lon(:), f_qfluxt_lat(:), & ! scalar flux reconstructed at cell centers
       f_epot(:), f_ekin(:), f_enthalpy(:), & ! time-averaged potential E, kinetic E and enthalpy
       f_epotfluxt(:), f_ekinfluxt(:), f_enthalpyfluxt(:) ! time averaged 'fluxes' of epot, ekin and enthalpy
  LOGICAL :: diagflux_on
  !$OMP THREADPRIVATE(diagflux_on)

  PUBLIC :: diagflux_on, init_diagflux, zero_qfluxt, qflux_centered_lonlat, flux_centered_lonlat, diagflux_energy

CONTAINS

  SUBROUTINE init_diagflux
    USE getin_mod
    INTEGER :: ll
    diagflux_on = .FALSE.
    CALL getin("diagflux", diagflux_on)
    ll = MERGE(llm,1,diagflux_on)
    CALL allocate_field(f_masst,         field_t,type_real,ll,       name="masst")
    CALL allocate_field(f_epot,          field_t,type_real,ll,       name="epot")
    CALL allocate_field(f_ekin,          field_t,type_real,ll,       name="ekin")
    CALL allocate_field(f_enthalpy,      field_t,type_real,ll,       name="enthalpy")
    CALL allocate_field(f_qmasst,        field_t,type_real,ll,nqtot, name="qmasst")
    CALL allocate_field(f_massfluxt,     field_u,type_real,ll,       name="massfluxt")
    CALL allocate_field(f_epotfluxt,     field_u,type_real,ll,       name="epotfluxt")
    CALL allocate_field(f_ekinfluxt,     field_u,type_real,ll,       name="ekinfluxt")
    CALL allocate_field(f_enthalpyfluxt, field_u,type_real,ll,       name="enthalpyfluxt")
    CALL allocate_field(f_qfluxt,        field_u,type_real,ll,nqtot, name="qfluxt")
    CALL allocate_field(f_qfluxt_lon,    field_t,type_real,ll,nqtot, name="qfluxt_lon")
    CALL allocate_field(f_qfluxt_lat,    field_t,type_real,ll,nqtot, name="qfluxt_lat")
    IF(diagflux_on) CALL zero_qfluxt
  END SUBROUTINE init_diagflux

#define ZERO2(field) buf2=field(ind) ; buf2(:,ll_begin:ll_end)=0.
#define ZERO3(field) buf3=field(ind) ; buf3(:,ll_begin:ll_end,:)=0.

  SUBROUTINE zero_qfluxt
    INTEGER :: ind
    REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       ZERO2(f_masst)
       ZERO2(f_epot)
       ZERO2(f_ekin)
       ZERO2(f_enthalpy)
       ZERO3(f_qmasst)
       ZERO2(f_massfluxt)
       ZERO2(f_epotfluxt)
       ZERO2(f_ekinfluxt)
       ZERO2(f_enthalpyfluxt)
       ZERO3(f_qfluxt)
    END DO
  END SUBROUTINE zero_qfluxt

!------------------------------------ Reconstruct fluxes at cell centers ---------------------------------------

  SUBROUTINE qflux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
    REAL(rstd), INTENT(IN) :: scale
    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:)
    INTEGER :: ind, itrac
    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       flux=f_flux(ind)
       flux_lon=f_flux_lon(ind)
       flux_lat=f_flux_lat(ind)
       DO itrac=1,nqtot
          CALL compute_flux_centered_lonlat(scale, flux(:,:,itrac), flux_lon(:,:,itrac), flux_lat(:,:,itrac))
       END DO
    END DO
  END SUBROUTINE qflux_centered_lonlat
  
  SUBROUTINE flux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
    REAL(rstd), INTENT(IN) :: scale
    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
    REAL(rstd), POINTER :: flux(:,:), flux_lon(:,:), flux_lat(:,:)
    INTEGER :: ind
    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       flux=f_flux(ind)
       flux_lon=f_flux_lon(ind)
       flux_lat=f_flux_lat(ind)
       CALL compute_flux_centered_lonlat(scale, flux, flux_lon, flux_lat)
    END DO
  END SUBROUTINE flux_centered_lonlat
  
  SUBROUTINE compute_flux_centered_lonlat(scale, flux, flux_lon, flux_lat)
    USE wind_mod
    REAL(rstd), INTENT(IN) :: scale
    REAL(rstd), INTENT(IN) :: flux(3*iim*jjm,llm)
    REAL(rstd), INTENT(OUT) :: flux_lon(iim*jjm,llm), flux_lat(iim*jjm,llm)
    REAL(rstd) :: flux_3d(iim*jjm,llm,3)
    CALL compute_flux_centered(scale, flux, flux_3d)
    CALL compute_wind_centered_lonlat_compound(flux_3d, flux_lon, flux_lat)
  END SUBROUTINE compute_flux_centered_lonlat

!------------------------------------ Compute energy fluxes ---------------------------------------

  SUBROUTINE diagflux_energy(frac, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
    REAL(rstd), INTENT(IN) :: frac
    TYPE(t_field),POINTER :: f_phis(:),f_rhodz(:),f_theta_rhodz(:),f_u(:), f_geopot(:), f_theta(:), f_hfluxt(:)
    REAL(rstd), POINTER :: phis(:), rhodz(:,:), theta_rhodz(:,:,:), u(:,:), &
         geopot(:,:), pk(:,:,:), hfluxt(:,:), &
         epot(:,:), ekin(:,:), enthalpy(:,:), &
         epotflux(:,:), ekinflux(:,:), enthalpyflux(:,:)
    INTEGER :: ind
    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       hfluxt = f_hfluxt(ind)
       phis = f_phis(ind)
       rhodz = f_rhodz(ind)
       theta_rhodz = f_theta_rhodz(ind)
       u = f_u(ind)
       geopot = f_geopot(ind)
       pk = f_theta(ind) ! buffer
       epot = f_epot(ind)
       ekin = f_ekin(ind)
       enthalpy = f_enthalpy(ind)
       epotflux = f_epotfluxt(ind)
       ekinflux = f_ekinfluxt(ind)
       enthalpyflux = f_enthalpyfluxt(ind)
       CALL compute_diagflux_energy(frac,hfluxt, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epotflux, ekinflux, enthalpyflux)
    END DO
  END SUBROUTINE diagflux_energy

  SUBROUTINE compute_diagflux_energy(frac, massflux, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epot_flux, ekin_flux, enthalpy_flux)
    USE disvert_mod, ONLY : ptop
    REAL(rstd), INTENT(IN) :: frac
    REAL(rstd), INTENT(IN) :: massflux(3*iim*jjm,llm), u(3*iim*jjm,llm),&
                              phis(iim*jjm), rhodz(iim*jjm,llm), theta_rhodz(iim*jjm,llm,nqtot)
    REAL(rstd), INTENT(INOUT) :: geopot(iim*jjm,llm+1), pk(iim*jjm,llm) ! pk = buffer
    REAL(rstd), INTENT(INOUT), DIMENSION(iim*jjm, llm)   ::  epot, ekin, enthalpy
    REAL(rstd), INTENT(INOUT), DIMENSION(3*iim*jjm, llm) ::  epot_flux, ekin_flux, enthalpy_flux    
    REAL(rstd) :: energy, p_ik, theta_ik, temp_ik, gv, Rd 
    INTEGER :: ij, l, ij_omp_begin_ext, ij_omp_end_ext
    Rd = kappa*cpp
    ! even if loops are of the _ext variant, we still need halo exchanges before reconstructing fluxes at cell centers
    ! => loop over interior region
    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin_ext,ij_omp_end_ext)
    ij_omp_begin_ext = ij_omp_begin_ext+ij_begin-1
    ij_omp_end_ext = ij_omp_end_ext+ij_begin-1
#include "../kernels/energy_fluxes.k90"
  END SUBROUTINE compute_diagflux_energy

END MODULE diagflux_mod