source: codes/icosagcm/trunk/src/kernels/energy_fluxes.k90 @ 599

Last change on this file since 599 was 599, checked in by dubos, 7 years ago

trunk : backported commits r582-r598 (transport diagnostics)
File size: 4.1 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- energy_fluxes ----------------------------------
3   ! First diagnose geopotential and temperature, column-wise
4   !$OMP BARRIER
5   DO ij=ij_omp_begin_ext,ij_omp_end_ext
6      pk(ij,llm) = ptop + .5*g*rhodz(ij,llm)
7   END DO
8   DO l = llm-1,1,-1
9      DO ij=ij_omp_begin_ext,ij_omp_end_ext
10         pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l)+rhodz(ij,l+1) )
11      END DO
12   END DO
13   SELECT CASE(caldyn_thermo)
14   CASE(thermo_theta)
15      DO l = 1,llm
16         DO ij=ij_omp_begin_ext,ij_omp_end_ext
17            p_ik = pk(ij,l)
18            theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
19            temp_ik = theta_ik*(p_ik/preff)**kappa
20            gv = (g*Rd)*temp_ik/p_ik
21            pk(ij,l) = temp_ik
22            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
23         END DO
24      END DO
25   CASE(thermo_entropy)
26      DO l = 1,llm
27         DO ij=ij_omp_begin_ext,ij_omp_end_ext
28            p_ik = pk(ij,l)
29            theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
30            temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
31            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
32            pk(ij,l) = temp_ik
33            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
34         END DO
35      END DO
36   END SELECT
37   !$OMP BARRIER
38   ! Now accumulate energies and energy fluxes
39   ! enthalpy
40   ! NB : at this point temperature is stored in array pk
41   ! pk then serves as buffer to store energy
42   DO l = ll_begin, ll_end
43      !DIR$ SIMD
44      DO ij=ij_begin_ext, ij_end_ext
45         energy = cpp*pk(ij,l)
46         enthalpy(ij,l) = enthalpy(ij,l) + frac*rhodz(ij,l)*energy
47         pk(ij,l) = energy
48      END DO
49   END DO
50   DO l = ll_begin, ll_end
51      !DIR$ SIMD
52      DO ij=ij_begin_ext, ij_end_ext
53         enthalpy_flux(ij+u_right,l) = enthalpy_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
54         enthalpy_flux(ij+u_lup,l) = enthalpy_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
55         enthalpy_flux(ij+u_ldown,l) = enthalpy_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
56      END DO
57   END DO
58   ! potential energy
59   DO l = ll_begin, ll_end
60      !DIR$ SIMD
61      DO ij=ij_begin_ext, ij_end_ext
62         energy = .5*(geopot(ij,l+1)+geopot(ij,l))
63         epot(ij,l) = epot(ij,l) + frac*rhodz(ij,l)*energy
64         pk(ij,l) = energy
65      END DO
66   END DO
67   DO l = ll_begin, ll_end
68      !DIR$ SIMD
69      DO ij=ij_begin_ext, ij_end_ext
70         epot_flux(ij+u_right,l) = epot_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
71         epot_flux(ij+u_lup,l) = epot_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
72         epot_flux(ij+u_ldown,l) = epot_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
73      END DO
74   END DO
75   ! kinetic energy
76   DO l = ll_begin, ll_end
77      !DIR$ SIMD
78      DO ij=ij_begin_ext, ij_end_ext
79         energy=0.d0
80         energy = energy + le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2
81         energy = energy + le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2
82         energy = energy + le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2
83         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2
84         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2
85         energy = energy + le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2
86         energy = energy * (.25/Ai(ij))
87         ekin(ij,l) = ekin(ij,l) + frac*rhodz(ij,l)*energy
88         pk(ij,l) = energy
89      END DO
90   END DO
91   DO l = ll_begin, ll_end
92      !DIR$ SIMD
93      DO ij=ij_begin_ext, ij_end_ext
94         ekin_flux(ij+u_right,l) = ekin_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
95         ekin_flux(ij+u_lup,l) = ekin_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
96         ekin_flux(ij+u_ldown,l) = ekin_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
97      END DO
98   END DO
99   !---------------------------- energy_fluxes ----------------------------------
100   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.