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

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

trunk : backported r600-603 from devel
File size: 7.4 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   ! NB : at this point pressure is stored in array pk
14   ! pk then serves as buffer to store temperature
15   SELECT CASE(caldyn_thermo)
16   CASE(thermo_theta)
17      DO l = 1,llm
18         DO ij=ij_omp_begin_ext,ij_omp_end_ext
19            p_ik = pk(ij,l)
20            theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
21            theta(ij,l) = theta_ik
22            temp_ik = theta_ik*(p_ik/preff)**kappa
23            gv = (g*Rd)*temp_ik/p_ik
24            pk(ij,l) = temp_ik
25            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
26         END DO
27      END DO
28   CASE(thermo_entropy)
29      DO l = 1,llm
30         DO ij=ij_omp_begin_ext,ij_omp_end_ext
31            p_ik = pk(ij,l)
32            theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
33            temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
34            theta(ij,l) = Treff*exp(theta_ik/cpp)
35            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
36            pk(ij,l) = temp_ik
37            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
38         END DO
39      END DO
40   END SELECT
41   !$OMP BARRIER
42   ! Now accumulate energies and energy fluxes
43   ! NB : at this point temperature is stored in array pk
44   ! pk then serves as buffer to store energy
45   ! enthalpy
46   DO l = ll_begin, ll_end
47      !DIR$ SIMD
48      DO ij=ij_begin_ext, ij_end_ext
49         energy = cpp*pk(ij,l)
50         enthalpy(ij,l) = enthalpy(ij,l) + frac*rhodz(ij,l)*energy
51         pk(ij,l) = energy
52      END DO
53   END DO
54   DO l = ll_begin, ll_end
55      !DIR$ SIMD
56      DO ij=ij_begin_ext, ij_end_ext
57         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))
58         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))
59         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))
60      END DO
61   END DO
62   ! potential energy
63   DO l = ll_begin, ll_end
64      !DIR$ SIMD
65      DO ij=ij_begin_ext, ij_end_ext
66         energy = .5*(geopot(ij,l+1)+geopot(ij,l))
67         epot(ij,l) = epot(ij,l) + frac*rhodz(ij,l)*energy
68         pk(ij,l) = energy
69      END DO
70   END DO
71   DO l = ll_begin, ll_end
72      !DIR$ SIMD
73      DO ij=ij_begin_ext, ij_end_ext
74         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))
75         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))
76         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))
77      END DO
78   END DO
79   ! theta
80   DO l = ll_begin, ll_end
81      !DIR$ SIMD
82      DO ij=ij_begin_ext, ij_end_ext
83         energy = theta(ij,l)
84         thetat(ij,l) = thetat(ij,l) + frac*rhodz(ij,l)*energy
85         pk(ij,l) = energy
86      END DO
87   END DO
88   DO l = ll_begin, ll_end
89      !DIR$ SIMD
90      DO ij=ij_begin_ext, ij_end_ext
91         thetat_flux(ij+u_right,l) = thetat_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
92         thetat_flux(ij+u_lup,l) = thetat_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
93         thetat_flux(ij+u_ldown,l) = thetat_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
94      END DO
95   END DO
96   ! kinetic energy
97   DO l = ll_begin, ll_end
98      !DIR$ SIMD
99      DO ij=ij_begin_ext, ij_end_ext
100         energy=0.d0
101         energy = energy + le(ij+u_rup)*de(ij+u_rup)*ue(ij+u_rup,l)**2
102         energy = energy + le(ij+u_lup)*de(ij+u_lup)*ue(ij+u_lup,l)**2
103         energy = energy + le(ij+u_left)*de(ij+u_left)*ue(ij+u_left,l)**2
104         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*ue(ij+u_ldown,l)**2
105         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*ue(ij+u_rdown,l)**2
106         energy = energy + le(ij+u_right)*de(ij+u_right)*ue(ij+u_right,l)**2
107         energy = energy * (.25/Ai(ij))
108         ekin(ij,l) = ekin(ij,l) + frac*rhodz(ij,l)*energy
109         pk(ij,l) = energy
110      END DO
111   END DO
112   DO l = ll_begin, ll_end
113      !DIR$ SIMD
114      DO ij=ij_begin_ext, ij_end_ext
115         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))
116         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))
117         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))
118      END DO
119   END DO
120   ! ulon
121   DO l = ll_begin, ll_end
122      !DIR$ SIMD
123      DO ij=ij_begin_ext, ij_end_ext
124         cx=centroid(ij,1)
125         cy=centroid(ij,2)
126         cz=centroid(ij,3)
127         ux=0. ; uy=0. ; uz=0.
128         ue_le = ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
129         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_up,1))-cx)
130         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_up,2))-cy)
131         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_up,3))-cz)
132         ue_le = ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
133         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_up,1))-cx)
134         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_up,2))-cy)
135         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_up,3))-cz)
136         ue_le = ne_left*ue(ij+u_left,l)*le(ij+u_left)
137         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_ldown,1))-cx)
138         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_ldown,2))-cy)
139         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_ldown,3))-cz)
140         ue_le = ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
141         ux = ux + ue_le*(.5*(xyz_v(ij+z_ldown,1)+xyz_v(ij+z_down,1))-cx)
142         uy = uy + ue_le*(.5*(xyz_v(ij+z_ldown,2)+xyz_v(ij+z_down,2))-cy)
143         uz = uz + ue_le*(.5*(xyz_v(ij+z_ldown,3)+xyz_v(ij+z_down,3))-cz)
144         ue_le = ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
145         ux = ux + ue_le*(.5*(xyz_v(ij+z_rdown,1)+xyz_v(ij+z_down,1))-cx)
146         uy = uy + ue_le*(.5*(xyz_v(ij+z_rdown,2)+xyz_v(ij+z_down,2))-cy)
147         uz = uz + ue_le*(.5*(xyz_v(ij+z_rdown,3)+xyz_v(ij+z_down,3))-cz)
148         ue_le = ne_right*ue(ij+u_right,l)*le(ij+u_right)
149         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_rdown,1))-cx)
150         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_rdown,2))-cy)
151         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_rdown,3))-cz)
152         ulon_i = ux*elon_i(ij,1) + uy*elon_i(ij,2) + uz*elon_i(ij,3)
153         energy = ulon_i*(1./Ai(ij))
154         ulon(ij,l) = ulon(ij,l) + frac*rhodz(ij,l)*energy
155         pk(ij,l) = energy
156      END DO
157   END DO
158   DO l = ll_begin, ll_end
159      !DIR$ SIMD
160      DO ij=ij_begin_ext, ij_end_ext
161         ulon_flux(ij+u_right,l) = ulon_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
162         ulon_flux(ij+u_lup,l) = ulon_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
163         ulon_flux(ij+u_ldown,l) = ulon_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
164      END DO
165   END DO
166   !---------------------------- energy_fluxes ----------------------------------
167   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.