source: LMDZ6/branches/Amaury_dev/libf/phylmd/ener_conserv.F90

Last change on this file was 5158, checked in by abarral, 7 weeks ago

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