source: LMDZ6/trunk/libf/phylmd/ener_conserv.f90 @ 5875

Last change on this file since 5875 was 5819, checked in by rkazeroni, 2 months ago

For GPU porting of ener_conserv routine:

  • Put routine into module (speeds up source-to-source transformation)
  • Add "horizontal" comment to specify possible names of horizontal variables
  • Modernize declarations of characters for GPU computing
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 9.3 KB
Line 
1!$gpum horizontal klon
2MODULE ener_conserv_mod
3
4  PRIVATE
5
6  PUBLIC ener_conserv
7
8  CONTAINS
9
10subroutine ener_conserv(klon,klev,pdtphys, &
11 &                      puo,pvo,pto,qx,ivap,iliq,isol, &
12 &                      pun,pvn,ptn,pqn,pqln,pqsn,dtke,masse,exner,d_t_ec)
13
14!=============================================================
15! Energy conservation
16! Based on the TKE equation
17! The M2 and N2 terms at the origin of TKE production are
18! concerted into heating in the d_t_ec term
19! Option 1 is the standard
20!        101 is for M2 term only
21!        101 for N2 term only
22!         -1 is a previours treatment for kinetic energy only
23!  FH (hourdin@lmd.jussieu.fr), 2013/04/25
24!=============================================================
25
26!=============================================================
27! Declarations
28!=============================================================
29
30! From module
31USE compbl_mod_h
32USE yoethf_mod_h
33USE clesphys_mod_h
34USE phys_local_var_mod, ONLY : d_u_vdf,d_v_vdf,d_t_vdf,d_u_ajs,d_v_ajs,d_t_ajs, &
35 &                             d_u_con,d_v_con,d_t_con,d_t_diss
36USE phys_local_var_mod, ONLY : d_t_eva,d_t_lsc,d_q_eva,d_q_lsc
37USE phys_local_var_mod, ONLY : d_u_oro,d_v_oro,d_u_lif,d_v_lif
38USE phys_local_var_mod, ONLY : du_gwd_hines,dv_gwd_hines,dv_gwd_front,dv_gwd_rando
39USE phys_state_var_mod, ONLY : du_gwd_front,du_gwd_rando
40USE phys_output_var_mod, ONLY : bils_ec,bils_ech,bils_tke,bils_kinetic,bils_enthalp,bils_latent,bils_diss
41USE add_phys_tend_mod, ONLY : fl_cor_ebil
42USE infotrac_phy, ONLY: nqtot
43
44
45USE yomcst_mod_h
46IMPLICIT none
47
48
49! Arguments
50INTEGER, INTENT(IN) :: klon,klev
51REAL, INTENT(IN) :: pdtphys
52REAL, DIMENSION(klon,klev), INTENT(IN)      :: puo,pvo,pto
53REAL, DIMENSION(klon,klev,nqtot), INTENT(IN):: qx
54INTEGER, INTENT(IN)                         :: ivap, iliq, isol
55REAL, DIMENSION(klon,klev), INTENT(IN)      :: pun,pvn,ptn,pqn,pqln,pqsn
56REAL, DIMENSION(klon,klev), INTENT(IN)      :: masse,exner
57REAL, DIMENSION(klon,klev+1), INTENT(IN)    :: dtke
58!
59REAL, DIMENSION(klon,klev), INTENT(OUT)     :: d_t_ec
60
61! Local
62      integer k,i
63REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt
64REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt
65REAL, DIMENSION(klon,klev) :: d_u,d_v,d_t,zv,zu,d_t_ech, pqo, pql0, pqs0
66REAL ZRCPD
67
68CHARACTER (LEN=80) :: abort_message
69CHARACTER (LEN=20), PARAMETER :: modname = 'ener_conser'
70
71d_t_ec(:,:)=0.
72
73IF(ivap == 0) CALL abort_physic (modname,'can''t run without water vapour',1)
74IF(iliq == 0) CALL abort_physic (modname,'can''t run without liquid water',1)
75pqo  = qx(:,:,ivap)
76pql0 = qx(:,:,iliq)
77IF(isol /= 0) pqs0 = qx(:,:,isol)
78
79IF (iflag_ener_conserv==-1) THEN
80!+jld ec_conser
81   DO k = 1, klev
82   DO i = 1, klon
83     IF (fl_cor_ebil .GT. 0) then
84       ZRCPD = RCPD*(1.0+RVTMP2*(pqn(i,k)+pqln(i,k)+pqsn(i,k)))
85     ELSE
86       ZRCPD = RCPD*(1.0+RVTMP2*pqn(i,k))
87     ENDIF
88     d_t_ec(i,k)=0.5/ZRCPD &
89 &     *(puo(i,k)**2+pvo(i,k)**2-pun(i,k)**2-pvn(i,k)**2)
90   ENDDO
91   ENDDO
92!-jld ec_conser
93
94
95
96ELSEIF (iflag_ener_conserv>=1) THEN
97
98   IF (iflag_ener_conserv<=2) THEN
99!     print*,'ener_conserv pbl=',iflag_pbl
100      IF (iflag_pbl>=20 .AND. iflag_pbl<=27) THEN !d_t_diss accounts for conserv
101         d_t(:,:)=d_t_ajs(:,:)   ! d_t_ajs = adjust + thermals
102         d_u(:,:)=d_u_ajs(:,:)+d_u_con(:,:)
103         d_v(:,:)=d_v_ajs(:,:)+d_v_con(:,:)
104      ELSE
105         d_t(:,:)=d_t_vdf(:,:)+d_t_ajs(:,:)   ! d_t_ajs = adjust + thermals
106         d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)
107         d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)
108      ENDIF
109   ELSEIF (iflag_ener_conserv==101) THEN
110      d_t(:,:)=0.
111      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)
112      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)
113   ELSEIF (iflag_ener_conserv==110) THEN
114      d_t(:,:)=d_t_vdf(:,:)+d_t_ajs(:,:)
115      d_u(:,:)=0.
116      d_v(:,:)=0.
117
118   ELSEIF (iflag_ener_conserv==3) THEN
119      d_t(:,:)=0.
120      d_u(:,:)=0.
121      d_v(:,:)=0.
122   ELSEIF (iflag_ener_conserv==4) THEN
123      d_t(:,:)=0.
124      d_u(:,:)=d_u_vdf(:,:)
125      d_v(:,:)=d_v_vdf(:,:)
126   ELSEIF (iflag_ener_conserv==5) THEN
127      d_t(:,:)=d_t_vdf(:,:)
128      d_u(:,:)=d_u_vdf(:,:)
129      d_v(:,:)=d_v_vdf(:,:)
130   ELSEIF (iflag_ener_conserv==6) THEN
131      d_t(:,:)=d_t_vdf(:,:)
132      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)
133      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)
134   ELSEIF (iflag_ener_conserv==7) THEN
135      d_t(:,:)=d_t_vdf(:,:)+d_t_ajs(:,:)
136      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)
137      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)
138   ELSEIF (iflag_ener_conserv==8) THEN
139      d_t(:,:)=d_t_vdf(:,:)
140      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)
141      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)
142   ELSEIF (iflag_ener_conserv==9) THEN
143      d_t(:,:)=d_t_vdf(:,:)
144      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)+d_u_oro(:,:)
145      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)+d_v_oro(:,:)
146   ELSEIF (iflag_ener_conserv==10) THEN
147      d_t(:,:)=d_t_vdf(:,:)
148      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)+d_u_oro(:,:)+d_u_lif(:,:)
149      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)+d_v_oro(:,:)+d_v_lif(:,:)
150   ELSEIF (iflag_ener_conserv==11) THEN
151      d_t(:,:)=d_t_vdf(:,:)
152      d_u(:,:)=d_u_vdf(:,:)+d_u_ajs(:,:)+d_u_con(:,:)+d_u_oro(:,:)+d_u_lif(:,:)
153      d_v(:,:)=d_v_vdf(:,:)+d_v_ajs(:,:)+d_v_con(:,:)+d_v_oro(:,:)+d_v_lif(:,:)
154      IF (ok_hines) THEN
155         d_u_vdf(:,:)=d_u_vdf(:,:)+du_gwd_hines(:,:)
156         d_v_vdf(:,:)=d_v_vdf(:,:)+dv_gwd_hines(:,:)
157      ENDIF
158      IF (.not. ok_hines .and. ok_gwd_rando) THEN
159         d_u_vdf(:,:)=d_u_vdf(:,:)+du_gwd_front(:,:)
160         d_v_vdf(:,:)=d_v_vdf(:,:)+dv_gwd_front(:,:)
161      ENDIF
162      IF (ok_gwd_rando) THEN
163         d_u_vdf(:,:)=d_u_vdf(:,:)+du_gwd_rando(:,:)
164         d_v_vdf(:,:)=d_v_vdf(:,:)+dv_gwd_rando(:,:)
165      ENDIF
166   ELSE
167      abort_message = 'iflag_ener_conserv non prevu'
168      CALL abort_physic (modname,abort_message,1)
169   ENDIF
170
171!----------------------------------------------------------------------------
172! Two options wether we consider time integration in the energy conservation
173!----------------------------------------------------------------------------
174
175   if (iflag_ener_conserv==2) then
176      zu(:,:)=puo(:,:)
177      zv(:,:)=pvo(:,:)
178   else
179      IF (iflag_pbl>=20 .AND. iflag_pbl<=27) THEN
180         zu(:,:)=puo(:,:)+d_u_vdf(:,:)+0.5*d_u(:,:)
181         zv(:,:)=pvo(:,:)+d_v_vdf(:,:)+0.5*d_v(:,:)
182      ELSE
183         zu(:,:)=puo(:,:)+0.5*d_u(:,:)
184         zv(:,:)=pvo(:,:)+0.5*d_v(:,:)
185      ENDIF
186   endif
187
188   fluxu(:,klev+1)=0.
189   fluxv(:,klev+1)=0.
190   fluxt(:,klev+1)=0.
191
192   do k=klev,1,-1
193      fluxu(:,k)=fluxu(:,k+1)+masse(:,k)*d_u(:,k)
194      fluxv(:,k)=fluxv(:,k+1)+masse(:,k)*d_v(:,k)
195      fluxt(:,k)=fluxt(:,k+1)+masse(:,k)*d_t(:,k)/exner(:,k)
196   enddo
197
198   dddu(:,1)=2*zu(:,1)*fluxu(:,1)
199   dddv(:,1)=2*zv(:,1)*fluxv(:,1)
200   dddt(:,1)=(exner(:,1)-1.)*fluxt(:,1)
201
202   do k=2,klev
203      dddu(:,k)=(zu(:,k)-zu(:,k-1))*fluxu(:,k)
204      dddv(:,k)=(zv(:,k)-zv(:,k-1))*fluxv(:,k)
205      dddt(:,k)=(exner(:,k)-exner(:,k-1))*fluxt(:,k)
206   enddo
207   dddu(:,klev+1)=0.
208   dddv(:,klev+1)=0.
209   dddt(:,klev+1)=0.
210
211   do k=1,klev
212      d_t_ech(:,k)=-(rcpd*(dddt(:,k)+dddt(:,k+1)))/(2.*rcpd*masse(:,k))
213      d_t_ec(:,k)=-(dddu(:,k)+dddu(:,k+1)+dddv(:,k)+dddv(:,k+1))/(2.*rcpd*masse(:,k))+d_t_ech(:,k)
214   enddo
215
216ENDIF
217
218!================================================================
219!  Computation of integrated enthalpie and kinetic energy variation
220!  FH (hourdin@lmd.jussieu.fr), 2013/04/25
221!  bils_ec : energie conservation term
222!  bils_ech : part of this term linked to temperature
223!  bils_tke : change of TKE
224!  bils_diss : dissipation of TKE (when activated)
225!  bils_kinetic : change of kinetic energie of the column
226!  bils_enthalp : change of enthalpie
227!  bils_latent  : change of latent heat. Computed between
228!          after reevaporation (at the beginning of the physics)
229!          and before large scale condensation (fisrtilp)
230!================================================================
231
232      bils_ec(:)=0.
233      bils_ech(:)=0.
234      bils_tke(:)=0.
235      bils_diss(:)=0.
236      bils_kinetic(:)=0.
237      bils_enthalp(:)=0.
238      bils_latent(:)=0.
239      DO k=1,klev
240        bils_ec(:)=bils_ec(:)-d_t_ec(:,k)*masse(:,k)
241        bils_diss(:)=bils_diss(:)-d_t_diss(:,k)*masse(:,k)
242        bils_kinetic(:)=bils_kinetic(:)+masse(:,k)* &
243     &           (pun(:,k)*pun(:,k)+pvn(:,k)*pvn(:,k) &
244     &            -puo(:,k)*puo(:,k)-pvo(:,k)*pvo(:,k))
245        bils_enthalp(:)= &
246     &  bils_enthalp(:)+masse(:,k)*(ptn(:,k)-pto(:,k)+d_t_ec(:,k)-d_t_eva(:,k)-d_t_lsc(:,k))
247!    &  bils_enthalp(:)+masse(:,k)*(ptn(:,k)-pto(:,k)+d_t_ec(:,k))
248        bils_latent(:)=bils_latent(:)+masse(:,k)* &
249!    &             (pqn(:,k)-pqo(:,k))
250     &             (pqn(:,k)-pqo(:,k)-d_q_eva(:,k)-d_q_lsc(:,k))
251      ENDDO
252      bils_ec(:)=rcpd*bils_ec(:)/pdtphys
253      bils_diss(:)=rcpd*bils_diss(:)/pdtphys
254      bils_kinetic(:)= 0.5*bils_kinetic(:)/pdtphys
255      bils_enthalp(:)=rcpd*bils_enthalp(:)/pdtphys
256      bils_latent(:)=rlvtt*bils_latent(:)/pdtphys
257!jyg<
258      IF (iflag_pbl > 1) THEN
259        DO k=1,klev
260          bils_tke(:)=bils_tke(:)+0.5*(dtke(:,k)+dtke(:,k+1))*masse(:,k)
261        ENDDO
262        bils_tke(:)=bils_tke(:)/pdtphys
263      ENDIF  ! (iflag_pbl > 1)
264!>jyg
265
266IF (iflag_ener_conserv>=1) THEN
267      bils_ech(:)=0.
268      DO k=1,klev
269        bils_ech(:)=bils_ech(:)-d_t_ech(:,k)*masse(:,k)
270      ENDDO
271      bils_ech(:)=rcpd*bils_ech(:)/pdtphys
272ENDIF
273
274RETURN
275
276END SUBROUTINE ener_conserv
277
278END MODULE ener_conserv_mod
Note: See TracBrowser for help on using the repository browser.