source: codes/icosagcm/trunk/src/diagnostics/diagflux.F90

Last change on this file was 953, checked in by adurocher, 5 years ago

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )
File size: 7.9 KB
Line 
1MODULE diagflux_mod
2  USE icosa
3  USE omp_para
4  IMPLICIT NONE
5  SAVE
6  PRIVATE
7
8  TYPE(t_field), POINTER, PUBLIC :: &
9       f_masst(:), f_qmasst(:), & ! time-averaged mass, tracer mass,
10       f_massfluxt(:), f_qfluxt(:), & ! time-integrated mass flux and tracer flux
11       f_qfluxt_lon(:), f_qfluxt_lat(:), & ! scalar flux reconstructed at cell centers
12       f_ulont(:), f_thetat(:), f_epot(:), f_ekin(:), f_enthalpy(:), & ! time-averaged potential E, kinetic E and enthalpy
13       f_ulonfluxt(:), f_thetafluxt(:), f_epotfluxt(:), f_ekinfluxt(:), f_enthalpyfluxt(:) ! time averaged 'fluxes' of epot, ekin and enthalpy
14  LOGICAL :: diagflux_on
15  !$OMP THREADPRIVATE(diagflux_on)
16
17  PUBLIC :: diagflux_on, init_diagflux, zero_qfluxt, qflux_centered_lonlat, flux_centered_lonlat, diagflux_energy
18
19CONTAINS
20
21  SUBROUTINE init_diagflux
22    USE getin_mod
23    USE abort_mod
24    INTEGER :: ll
25    diagflux_on = .FALSE.
26    CALL getin("diagflux", diagflux_on)
27    IF (diagflux_on) THEN
28       CALL abort_acc("diagflux /= .FALSE.")
29    END IF
30
31    ll = MERGE(llm,1,diagflux_on)
32    CALL allocate_field(f_masst,         field_t,type_real,ll,       name="masst")
33    CALL allocate_field(f_ulont,         field_t,type_real,ll,       name="ulont")
34    CALL allocate_field(f_thetat,        field_t,type_real,ll,       name="thetat")
35    CALL allocate_field(f_epot,          field_t,type_real,ll,       name="epot")
36    CALL allocate_field(f_ekin,          field_t,type_real,ll,       name="ekin")
37    CALL allocate_field(f_enthalpy,      field_t,type_real,ll,       name="enthalpy")
38    CALL allocate_field(f_qmasst,        field_t,type_real,ll,nqtot, name="qmasst")
39    CALL allocate_field(f_massfluxt,     field_u,type_real,ll,       name="massfluxt")
40    CALL allocate_field(f_ulonfluxt,     field_u,type_real,ll,       name="ulonfluxt")
41    CALL allocate_field(f_thetafluxt,    field_u,type_real,ll,       name="thetafluxt")
42    CALL allocate_field(f_epotfluxt,     field_u,type_real,ll,       name="epotfluxt")
43    CALL allocate_field(f_ekinfluxt,     field_u,type_real,ll,       name="ekinfluxt")
44    CALL allocate_field(f_enthalpyfluxt, field_u,type_real,ll,       name="enthalpyfluxt")
45    CALL allocate_field(f_qfluxt,        field_u,type_real,ll,nqtot, name="qfluxt")
46    CALL allocate_field(f_qfluxt_lon,    field_t,type_real,ll,nqtot, name="qfluxt_lon")
47    CALL allocate_field(f_qfluxt_lat,    field_t,type_real,ll,nqtot, name="qfluxt_lat")
48    IF(diagflux_on) CALL zero_qfluxt
49  END SUBROUTINE init_diagflux
50
51#define ZERO2(field) buf2=field(ind) ; buf2(:,ll_begin:ll_end)=0.
52#define ZERO3(field) buf3=field(ind) ; buf3(:,ll_begin:ll_end,:)=0.
53
54  SUBROUTINE zero_qfluxt
55    INTEGER :: ind
56    REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
57    DO ind=1,ndomain
58       IF (.NOT. assigned_domain(ind)) CYCLE
59       CALL swap_dimensions(ind)
60       ZERO2(f_masst)
61       ZERO2(f_ulont)
62       ZERO2(f_thetat)
63       ZERO2(f_epot)
64       ZERO2(f_ekin)
65       ZERO2(f_enthalpy)
66       ZERO3(f_qmasst)
67       ZERO2(f_massfluxt)
68       ZERO2(f_ulonfluxt)
69       ZERO2(f_thetafluxt)
70       ZERO2(f_epotfluxt)
71       ZERO2(f_ekinfluxt)
72       ZERO2(f_enthalpyfluxt)
73       ZERO3(f_qfluxt)
74    END DO
75  END SUBROUTINE zero_qfluxt
76
77!------------------------------------ Reconstruct fluxes at cell centers ---------------------------------------
78
79  SUBROUTINE qflux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
80    REAL(rstd), INTENT(IN) :: scale
81    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
82    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:)
83    INTEGER :: ind, itrac
84    DO ind=1,ndomain
85       IF (.NOT. assigned_domain(ind)) CYCLE
86       CALL swap_dimensions(ind)
87       CALL swap_geometry(ind)
88       flux=f_flux(ind)
89       flux_lon=f_flux_lon(ind)
90       flux_lat=f_flux_lat(ind)
91       DO itrac=1,nqtot
92          CALL compute_flux_centered_lonlat(scale, flux(:,:,itrac), flux_lon(:,:,itrac), flux_lat(:,:,itrac))
93       END DO
94    END DO
95  END SUBROUTINE qflux_centered_lonlat
96 
97  SUBROUTINE flux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
98    REAL(rstd), INTENT(IN) :: scale
99    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
100    REAL(rstd), POINTER :: flux(:,:), flux_lon(:,:), flux_lat(:,:)
101    INTEGER :: ind
102    DO ind=1,ndomain
103       IF (.NOT. assigned_domain(ind)) CYCLE
104       CALL swap_dimensions(ind)
105       CALL swap_geometry(ind)
106       flux=f_flux(ind)
107       flux_lon=f_flux_lon(ind)
108       flux_lat=f_flux_lat(ind)
109       CALL compute_flux_centered_lonlat(scale, flux, flux_lon, flux_lat)
110    END DO
111  END SUBROUTINE flux_centered_lonlat
112 
113  SUBROUTINE compute_flux_centered_lonlat(scale, flux, flux_lon, flux_lat)
114    USE wind_mod
115    REAL(rstd), INTENT(IN) :: scale
116    REAL(rstd), INTENT(IN) :: flux(3*iim*jjm,llm)
117    REAL(rstd), INTENT(OUT) :: flux_lon(iim*jjm,llm), flux_lat(iim*jjm,llm)
118    REAL(rstd) :: flux_3d(iim*jjm,llm,3)
119    CALL compute_flux_centered(scale, flux, flux_3d)
120    CALL compute_wind_centered_lonlat_compound(flux_3d, flux_lon, flux_lat)
121  END SUBROUTINE compute_flux_centered_lonlat
122
123!------------------------------------ Compute energy fluxes ---------------------------------------
124
125  SUBROUTINE diagflux_energy(frac, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_pk, f_hfluxt)
126    REAL(rstd), INTENT(IN) :: frac
127    TYPE(t_field),POINTER :: f_phis(:),f_rhodz(:),f_theta_rhodz(:),f_u(:), f_geopot(:), f_theta(:), f_pk(:), f_hfluxt(:)
128    REAL(rstd), POINTER :: phis(:), rhodz(:,:), theta_rhodz(:,:,:), u(:,:), &
129         geopot(:,:), theta(:,:,:), pk(:,:), hfluxt(:,:), &
130         ulont(:,:), thetat(:,:), epot(:,:), ekin(:,:), enthalpy(:,:), &
131         thetaflux(:,:), ulonflux(:,:), epotflux(:,:), ekinflux(:,:), enthalpyflux(:,:)
132    INTEGER :: ind
133    DO ind=1,ndomain
134       IF (.NOT. assigned_domain(ind)) CYCLE
135       CALL swap_dimensions(ind)
136       CALL swap_geometry(ind)
137       hfluxt = f_hfluxt(ind)
138       phis = f_phis(ind)
139       rhodz = f_rhodz(ind)
140       theta_rhodz = f_theta_rhodz(ind)
141       u = f_u(ind)
142       geopot = f_geopot(ind)
143       theta = f_theta(ind) ! buffer
144       pk = f_pk(ind) ! buffer
145       ulont = f_ulont(ind)
146       thetat = f_thetat(ind)
147       epot = f_epot(ind)
148       ekin = f_ekin(ind)
149       enthalpy = f_enthalpy(ind)
150       ulonflux = f_ulonfluxt(ind)
151       thetaflux = f_thetafluxt(ind)
152       epotflux = f_epotfluxt(ind)
153       ekinflux = f_ekinfluxt(ind)
154       enthalpyflux = f_enthalpyfluxt(ind)
155       CALL compute_diagflux_energy(frac,hfluxt, phis,rhodz,theta_rhodz,u, geopot,theta,pk, &
156            ulont, thetat, epot, ekin, enthalpy, &
157            ulonflux, thetaflux, epotflux, ekinflux, enthalpyflux)
158    END DO
159  END SUBROUTINE diagflux_energy
160
161  SUBROUTINE compute_diagflux_energy(frac, massflux, phis,rhodz,theta_rhodz,ue, geopot,theta,pk, &
162       ulon, thetat, epot, ekin, enthalpy, &
163       ulon_flux, thetat_flux, epot_flux, ekin_flux, enthalpy_flux)
164    USE disvert_mod, ONLY : ptop
165    REAL(rstd), INTENT(IN) :: frac
166    REAL(rstd), INTENT(IN) :: massflux(3*iim*jjm,llm), ue(3*iim*jjm,llm),&
167                              phis(iim*jjm), rhodz(iim*jjm,llm), theta_rhodz(iim*jjm,llm,nqtot)
168    REAL(rstd), INTENT(INOUT) :: geopot(iim*jjm,llm+1), theta(iim*jjm,llm), pk(iim*jjm,llm) ! theta,pk = buffers
169    REAL(rstd), INTENT(INOUT), DIMENSION(iim*jjm, llm)   ::  ulon, thetat, epot, ekin, enthalpy
170    REAL(rstd), INTENT(INOUT), DIMENSION(3*iim*jjm, llm) ::  ulon_flux, thetat_flux, epot_flux, ekin_flux, enthalpy_flux   
171    REAL(rstd) :: energy, p_ik, theta_ik, temp_ik, gv, Rd, cx,cy,cz, ux,uy,uz, ue_le,ulon_i
172    INTEGER :: ij, l, ij_omp_begin_ext, ij_omp_end_ext
173    Rd = kappa*cpp
174    ! even if loops are of the _ext variant, we still need halo exchanges before reconstructing fluxes at cell centers
175    ! => loop over interior region
176    CALL distrib_level(ij_begin_ext, ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
177#include "../kernels/energy_fluxes.k90"
178  END SUBROUTINE compute_diagflux_energy
179
180END MODULE diagflux_mod
Note: See TracBrowser for help on using the repository browser.