source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/add_phys_tend_mod.F90 @ 3324

Last change on this file since 3324 was 2848, checked in by jyg, 8 years ago

Changes in module add_phys_tend_mod.F90: (i)
vertical integrations are performed in subroutine
integr_v; (ii) the new subroutine diag_phys_tend
makes it possible to check water and energy
conservation of a set of tendencies added to any
atmospheric state.
File size: 29.5 KB
Line 
1!
2MODULE add_phys_tend_mod
3
4  IMPLICIT NONE 
5  ! flag to compute diagnostics to check energy conservation. If  fl_ebil==0, no check
6  INTEGER, SAVE ::   fl_ebil
7!$OMP THREADPRIVATE(fl_ebil)
8  ! flag to include modifcations to ensure energy conservation. If  fl_cor_ebil==0, no corrections
9  ! Note that with time, all these modifications should be included by default
10  INTEGER, SAVE ::   fl_cor_ebil
11!$OMP THREADPRIVATE(fl_cor_ebil)
12
13CONTAINS
14
15SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap)
16  ! ======================================================================
17  ! Ajoute les tendances de couche limite, soit determinees par la
18  ! parametrisation
19  ! physique, soit forcees,  aux variables d etat de la dynamique t_seri,
20  ! q_seri ...
21  ! ======================================================================
22
23
24  ! ======================================================================
25  ! Declarations
26  ! ======================================================================
27
28  USE dimphy, ONLY: klon, klev
29!  USE dimphy
30  USE phys_local_var_mod
31  USE phys_state_var_mod
32  USE mod_grid_phy_lmdz, ONLY: nbp_lev
33  IMPLICIT NONE
34  REAL,SAVE,ALLOCATABLE :: hthturb_gcssold(:)
35  REAL,SAVE,ALLOCATABLE :: hqturb_gcssold(:)
36!$OMP THREADPRIVATE(hthturb_gcssold,hqturb_gcssold)
37  REAL,SAVE :: dtime_frcg
38  LOGICAL,SAVE :: turb_fcg_gcssold
39  LOGICAL,SAVE :: firstcall=.true.
40!$OMP THREADPRIVATE(firstcall,turb_fcg_gcssold,dtime_frcg)
41  INTEGER abortphy
42!  COMMON /turb_forcing/dtime_frcg, hthturb_gcssold, hqturb_gcssold, &
43!    turb_fcg_gcssold
44
45  ! Arguments :
46  ! ------------
47  REAL zdu(klon, klev), zdv(klon, klev)
48  REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev)
49  CHARACTER *(*) text
50  REAL paprs(klon,klev+1)
51  INTEGER flag_inhib_tend ! if flag_inhib_tend != 0, tendencies are not added
52  INTEGER itap ! time step number
53
54  ! Local :
55  ! --------
56  REAL zzdt(klon, klev), zzdq(klon, klev)
57  INTEGER i, k
58
59  IF (firstcall) THEN
60    ALLOCATE(hthturb_gcssold(nbp_lev))
61    ALLOCATE(hqturb_gcssold(nbp_lev))
62    firstcall=.false.
63  ENDIF
64
65  IF (turb_fcg_gcssold) THEN
66    DO k = 1, klev
67      DO i = 1, klon
68        zzdt(i, k) = hthturb_gcssold(k)*dtime_frcg
69        zzdq(i, k) = hqturb_gcssold(k)*dtime_frcg
70      END DO
71    END DO
72    PRINT *, ' add_pbl_tend, dtime_frcg ', dtime_frcg
73    PRINT *, ' add_pbl_tend, zzdt ', zzdt
74    PRINT *, ' add_pbl_tend, zzdq ', zzdq
75    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0)
76  ELSE
77    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0)
78  END IF
79
80
81  RETURN
82END SUBROUTINE add_pbl_tend
83!
84! $Id: add_phys_tend.F90 2611 2016-08-03 15:41:26Z jyg $
85!
86SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text, &
87                          abortphy,flag_inhib_tend, itap, diag_mode)
88!======================================================================
89! Ajoute les tendances des variables physiques aux variables
90! d'etat de la dynamique t_seri, q_seri ...
91! On en profite pour faire des tests sur les tendances en question.
92!======================================================================
93
94
95!======================================================================
96! Declarations
97!======================================================================
98
99USE dimphy, ONLY: klon, klev
100USE phys_state_var_mod, ONLY : dtime
101USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri
102USE phys_state_var_mod, ONLY: ftsol
103USE geometry_mod, ONLY: longitude_deg, latitude_deg
104USE print_control_mod, ONLY: prt_level
105USE cmp_seri_mod
106USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
107  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
108IMPLICIT none
109  include "YOMCST.h"
110  include "clesphys.h"
111
112! Arguments :
113!------------
114REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdu, zdv
115REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdt, zdql, zdqi
116REAL, DIMENSION(klon,klev+1),   INTENT(IN)    :: paprs
117CHARACTER*(*),                  INTENT(IN)    :: text
118INTEGER,                        INTENT(IN)    :: abortphy
119INTEGER,                        INTENT(IN)    :: flag_inhib_tend ! if not 0, tendencies are not added
120INTEGER,                        INTENT(IN)    :: itap            ! time step number
121INTEGER,                        INTENT(IN)    :: diag_mode       ! 0 -> normal effective mode
122                                                                 ! 1 -> only conservation stats are computed
123!
124REAL, DIMENSION(klon,klev),     INTENT(INOUT) :: zdq
125
126! Local :
127!--------
128REAL zt,zq
129REAL zq_int, zqp_int, zq_new
130
131REAL zqp(klev)
132
133! Save variables, used in diagnostic mode (diag_mode=1).
134REAL, DIMENSION(klon,klev)   :: sav_u_seri, sav_v_seri
135REAL, DIMENSION(klon,klev)   :: sav_ql_seri, sav_qs_seri, sav_q_seri
136REAL, DIMENSION(klon,klev)   :: sav_t_seri
137REAL, DIMENSION(klon,klev)   :: sav_zdq
138!
139INTEGER i, k,j, n
140INTEGER jadrs(klon*klev), jbad
141INTEGER jqadrs(klon*klev), jqbad
142INTEGER kadrs(klon*klev)
143INTEGER kqadrs(klon*klev)
144
145LOGICAL done(klon)
146
147integer debug_level
148logical, save :: first=.true.
149!$OMP THREADPRIVATE(first)
150!
151!======================================================================
152! Variables for energy conservation tests
153!======================================================================
154!
155
156! zh_col-------  total enthalpy of vertical air column
157! (air with watter vapour, liquid and solid) (J/m2)
158! zh_dair_col--- total enthalpy of dry air (J/m2)
159! zh_qw_col----  total enthalpy of watter vapour (J/m2)
160! zh_ql_col----  total enthalpy of liquid watter (J/m2)
161! zh_qs_col----  total enthalpy of solid watter  (J/m2)
162! zqw_col------  total mass of watter vapour (kg/m2)
163! zql_col------  total mass of liquid watter (kg/m2)
164! zqs_col------  total mass of solid watter (kg/m2)
165! zek_col------  total kinetic energy (kg/m2)
166!
167REAL zairm(klon, klev) ! layer air mass (kg/m2)
168REAL zqw_col(klon,2)
169REAL zql_col(klon,2)
170REAL zqs_col(klon,2)
171REAL zek_col(klon,2)
172REAL zh_dair_col(klon,2)
173REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2)
174REAL zh_col(klon,2)
175
176REAL zcpvap, zcwat, zcice
177
178!======================================================================
179! Initialisations
180
181     IF (prt_level >= 5) then
182        write (*,*) "In add_phys_tend, after ",text
183        call flush
184     end if
185
186     ! if flag_inhib_tend != 0, tendencies are not added
187     IF (flag_inhib_tend /= 0) then
188        ! If requiered, diagnostics are shown
189        IF (flag_inhib_tend > 0) then
190           ! print some diagnostics if xxx_seri have changed
191           call cmp_seri(flag_inhib_tend,text)
192        END IF
193        RETURN ! on n ajoute pas les tendance
194     END IF
195
196     IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele
197                              ! a deja plante.
198
199     debug_level=10
200     if (first) then
201        print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs
202        first=.false.
203     endif
204!
205!  print *,'add_phys_tend: paprs ',paprs
206! When in diagnostic mode, save initial values of out variables
207  IF (diag_mode == 1) THEN
208    sav_u_seri(:,:)  = u_seri(:,:)
209    sav_v_seri(:,:)  = v_seri(:,:)
210    sav_ql_seri(:,:) = ql_seri(:,:)
211    sav_qs_seri(:,:) = qs_seri(:,:)
212    sav_q_seri(:,:)  = q_seri(:,:)
213    sav_t_seri(:,:)  = t_seri(:,:)
214    sav_zdq(:,:)     = zdq(:,:)   
215  ENDIF ! (diag_mode == 1)
216!======================================================================
217! Diagnostics for energy conservation tests
218!======================================================================
219  DO k = 1, klev
220    ! layer air mass
221    zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg
222  END DO
223
224  if (fl_ebil .GT. 0) then
225    ! ------------------------------------------------
226    ! Compute vertical sum for each atmospheric column
227    ! ------------------------------------------------
228    n=1   ! begining of time step
229
230    zcpvap = rcpv
231    zcwat = rcw
232    zcice = rcs
233
234    CALL integr_v(klon, klev, zcpvap, &
235                  t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
236                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
237                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
238
239  end if ! end if (fl_ebil .GT. 0)
240
241!======================================================================
242! Ajout des tendances sur le vent et l'eau liquide
243!======================================================================
244
245     u_seri(:,:)=u_seri(:,:)+zdu(:,:)
246     v_seri(:,:)=v_seri(:,:)+zdv(:,:)
247     ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
248     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
249
250!======================================================================
251! On ajoute les tendances de la temperature et de la vapeur d'eau
252! en verifiant que ca ne part pas dans les choux
253!======================================================================
254
255      jbad=0
256      jqbad=0
257      DO k = 1, klev
258         DO i = 1, klon
259            zt=t_seri(i,k)+zdt(i,k)
260            zq=q_seri(i,k)+zdq(i,k)
261            IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then
262            jbad = jbad + 1
263            jadrs(jbad) = i
264            kadrs(jbad) = k
265            ENDIF
266            IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
267            jqbad = jqbad + 1
268            jqadrs(jqbad) = i
269            kqadrs(jqbad) = k
270            ENDIF
271            t_seri(i,k)=zt
272            q_seri(i,k)=zq
273         ENDDO
274      ENDDO
275
276!=====================================================================================
277! Impression et stop en cas de probleme important
278!=====================================================================================
279
280IF (jbad .GT. 0) THEN
281      DO j = 1, jbad
282         i=jadrs(j)
283         if(prt_level.ge.debug_level) THEN
284          print*,'PLANTAGE POUR LE POINT i lon lat =',&
285                 i,longitude_deg(i),latitude_deg(i),text
286          print*,'l    T     dT       Q     dQ    '
287          DO k = 1, klev
288             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
289          ENDDO
290          call print_debug_phys(i,debug_level,text)
291         endif
292      ENDDO
293ENDIF
294!
295!=====================================================================================
296! Impression, warning et correction en cas de probleme moins important
297!=====================================================================================
298IF (jqbad .GT. 0) THEN
299      done(:) = .false.                         !jyg
300      DO j = 1, jqbad
301        i=jqadrs(j)
302          if(prt_level.ge.debug_level) THEN
303           print*,'WARNING  : EAU POUR LE POINT i lon lat =',&
304                  i,longitude_deg(i),latitude_deg(i),text
305           print*,'l    T     dT       Q     dQ    '
306           DO k = 1, klev
307              write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
308           ENDDO
309          endif
310          IF (ok_conserv_q) THEN
311!jyg<20140228 Corrections pour conservation de l'eau
312            IF (.NOT.done(i)) THEN                  !jyg
313              DO k = 1, klev
314                zqp(k) = max(q_seri(i,k),1.e-15)
315              ENDDO
316              zq_int  = 0.
317              zqp_int = 0.
318              DO k = 1, klev
319                zq_int  = zq_int  + q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
320                zqp_int = zqp_int + zqp(k)     *(paprs(i,k)-paprs(i,k+1))/Rg
321              ENDDO
322              if(prt_level.ge.debug_level) THEN
323               print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
324                                    i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
325              endif
326              DO k = 1, klev
327                zq_new = zqp(k)*zq_int/zqp_int
328                zdq(i,k) = zdq(i,k) + zq_new - q_seri(i,k)
329                q_seri(i,k) = zq_new
330              ENDDO
331              done(i) = .true.
332            ENDIF !(.NOT.done(i))
333          ELSE
334!jyg>
335            DO k = 1, klev
336              zq=q_seri(i,k)+zdq(i,k)
337              if (zq.lt.1.e-15) then
338                 if (q_seri(i,k).lt.1.e-15) then
339                  if(prt_level.ge.debug_level) THEN
340                   print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
341                  endif
342                  q_seri(i,k)=1.e-15
343                  zdq(i,k)=(1.e-15-q_seri(i,k))
344                 endif
345              endif
346!              zq=q_seri(i,k)+zdq(i,k)
347!              if (zq.lt.1.e-15) then
348!                 zdq(i,k)=(1.e-15-q_seri(i,k))
349!              endif
350            ENDDO
351!jyg<20140228
352          ENDIF   !  (ok_conserv_q)
353!jyg>
354      ENDDO ! j = 1, jqbad
355ENDIF
356!
357
358!IM ajout memes tests pour reverifier les jbad, jqbad beg
359      jbad=0
360      jqbad=0
361      DO k = 1, klev
362         DO i = 1, klon
363            IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
364            jbad = jbad + 1
365            jadrs(jbad) = i
366!            if(prt_level.ge.debug_level) THEN
367!             print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
368!            endif
369            ENDIF
370            IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
371            jqbad = jqbad + 1
372            jqadrs(jqbad) = i
373            kqadrs(jqbad) = k
374!            if(prt_level.ge.debug_level) THEN
375!             print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
376!            endif
377            ENDIF
378         ENDDO
379      ENDDO
380IF (jbad .GT. 0) THEN
381      DO j = 1, jbad
382         i=jadrs(j)
383         k=kadrs(j)
384         if(prt_level.ge.debug_level) THEN
385          print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',&
386                 i,itap,longitude_deg(i),latitude_deg(i),text,jbad, &
387       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
388!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
389          print*,'l    T     dT       Q     dQ    '
390          DO k = 1, klev
391             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
392          ENDDO
393          call print_debug_phys(i,debug_level,text)
394         endif
395      ENDDO
396ENDIF
397!
398IF (jqbad .GT. 0) THEN
399      DO j = 1, jqbad
400         i=jqadrs(j)
401         k=kqadrs(j)
402         if(prt_level.ge.debug_level) THEN
403          print*,'WARNING  : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',&
404                 i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,&
405       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
406!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
407          print*,'l    T     dT       Q     dQ    '
408          DO k = 1, klev
409            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
410          ENDDO
411          call print_debug_phys(i,debug_level,text)
412         endif
413      ENDDO
414ENDIF
415
416!======================================================================
417! Contrôle des min/max pour arrêt du modèle
418! Si le modele est en mode abortphy, on retire les tendances qu'on
419! vient d'ajouter. Pas exactement parce qu'on ne tient pas compte des
420! seuils.
421!======================================================================
422
423      CALL hgardfou(t_seri,ftsol,text,abortphy)
424      IF (abortphy==1) THEN
425        Print*,'ERROR ABORT hgardfou dans ',text
426! JLD pourquoi on ne modifie pas de meme t_seri et q_seri ?
427        u_seri(:,:)=u_seri(:,:)-zdu(:,:)
428        v_seri(:,:)=v_seri(:,:)-zdv(:,:)
429        ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
430        qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
431      ENDIF
432
433!======================================================================
434! Diagnostics for energy conservation tests
435!======================================================================
436
437  if (fl_ebil .GT. 0) then
438 
439    ! ------------------------------------------------
440    ! Compute vertical sum for each atmospheric column
441    ! ------------------------------------------------
442    n=2   ! end of time step
443
444    CALL integr_v(klon, klev, zcpvap, &
445                  t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
446                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
447                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
448
449    ! ------------------------------------------------
450    ! Compute the changes by unit of time
451    ! ------------------------------------------------
452
453    d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/dtime
454    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/dtime
455    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/dtime
456    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
457
458    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/dtime
459
460    d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/dtime
461    d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/dtime
462    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/dtime
463    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/dtime
464
465    d_h_col = (zh_col(:,2)-zh_col(:,1))/dtime
466
467  end if ! end if (fl_ebil .GT. 0)
468!
469! When in diagnostic mode, restore "out" variables to initial values.
470  IF (diag_mode == 1) THEN
471    u_seri(:,:)  = sav_u_seri(:,:)
472    v_seri(:,:)  = sav_v_seri(:,:)
473    ql_seri(:,:) = sav_ql_seri(:,:)
474    qs_seri(:,:) = sav_qs_seri(:,:)
475    q_seri(:,:)  = sav_q_seri(:,:)
476    t_seri(:,:)  = sav_t_seri(:,:)
477    zdq(:,:)     = sav_zdq(:,:)   
478  ENDIF ! (mode == 1)
479
480  RETURN
481END SUBROUTINE add_phys_tend
482
483SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, &
484                          zdu,zdv,zdt,zdq,zdql,zdqs,paprs,text)
485!======================================================================
486! Ajoute les tendances des variables physiques aux variables
487! d'etat de la dynamique t_seri, q_seri ...
488! On en profite pour faire des tests sur les tendances en question.
489!======================================================================
490
491
492!======================================================================
493! Declarations
494!======================================================================
495
496USE phys_state_var_mod, ONLY : dtime, ftsol
497USE geometry_mod, ONLY: longitude_deg, latitude_deg
498USE print_control_mod, ONLY: prt_level
499USE cmp_seri_mod
500USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
501  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
502IMPLICIT none
503  include "YOMCST.h"
504  include "clesphys.h"
505
506! Arguments :
507!------------
508INTEGER, INTENT(IN)                           :: nlon, nlev
509REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: uu, vv
510REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: temp, qv, ql, qs
511REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdu, zdv
512REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdt, zdq, zdql, zdqs
513REAL, DIMENSION(nlon,nlev+1),   INTENT(IN)    :: paprs
514CHARACTER*(*),                  INTENT(IN)    :: text
515
516! Local :
517!--------
518REAL, DIMENSION(nlon,nlev)      :: uu_n, vv_n
519REAL, DIMENSION(nlon,nlev)      :: temp_n, qv_n, ql_n, qs_n
520
521
522!
523INTEGER k, n
524
525integer debug_level
526logical, save :: first=.true.
527!$OMP THREADPRIVATE(first)
528!
529!======================================================================
530! Variables for energy conservation tests
531!======================================================================
532!
533
534! zh_col-------  total enthalpy of vertical air column
535! (air with watter vapour, liquid and solid) (J/m2)
536! zh_dair_col--- total enthalpy of dry air (J/m2)
537! zh_qw_col----  total enthalpy of watter vapour (J/m2)
538! zh_ql_col----  total enthalpy of liquid watter (J/m2)
539! zh_qs_col----  total enthalpy of solid watter  (J/m2)
540! zqw_col------  total mass of watter vapour (kg/m2)
541! zql_col------  total mass of liquid watter (kg/m2)
542! zqs_col------  total mass of solid watter (kg/m2)
543! zek_col------  total kinetic energy (kg/m2)
544!
545REAL zairm(nlon, nlev) ! layer air mass (kg/m2)
546REAL zqw_col(nlon,2)
547REAL zql_col(nlon,2)
548REAL zqs_col(nlon,2)
549REAL zek_col(nlon,2)
550REAL zh_dair_col(nlon,2)
551REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2)
552REAL zh_col(nlon,2)
553
554!======================================================================
555! Initialisations
556
557     IF (prt_level >= 5) then
558        write (*,*) "In diag_phys_tend, after ",text
559        call flush
560     end if
561
562     debug_level=10
563     if (first) then
564        print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs
565        first=.false.
566     endif
567!
568!  print *,'add_phys_tend: paprs ',paprs
569!======================================================================
570! Diagnostics for energy conservation tests
571!======================================================================
572  DO k = 1, nlev
573    ! layer air mass
574    zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg
575  END DO
576
577  if (fl_ebil .GT. 0) then
578    ! ------------------------------------------------
579    ! Compute vertical sum for each atmospheric column
580    ! ------------------------------------------------
581    n=1   ! begining of time step
582
583    CALL integr_v(nlon, nlev, rcpv, &
584                  temp, qv, ql, qs, uu, vv, zairm, &
585                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
586                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
587
588  end if ! end if (fl_ebil .GT. 0)
589
590!======================================================================
591! Ajout des tendances sur le vent, la temperature et les diverses phases de l'eau
592!======================================================================
593
594     uu_n(:,:)=uu(:,:)+zdu(:,:)
595     vv_n(:,:)=vv(:,:)+zdv(:,:)
596     qv_n(:,:)=qv(:,:)+zdq(:,:)
597     ql_n(:,:)=ql(:,:)+zdql(:,:)
598     qs_n(:,:)=qs(:,:)+zdqs(:,:)
599     temp_n(:,:)=temp(:,:)+zdt(:,:)
600
601
602
603!======================================================================
604! Diagnostics for energy conservation tests
605!======================================================================
606
607  if (fl_ebil .GT. 0) then
608 
609    ! ------------------------------------------------
610    ! Compute vertical sum for each atmospheric column
611    ! ------------------------------------------------
612    n=2   ! end of time step
613
614    CALL integr_v(nlon, nlev, rcpv, &
615                  temp_n, qv_n, ql_n, qs_n, uu_n, vv_n, zairm, &
616                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
617                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
618
619    ! ------------------------------------------------
620    ! Compute the changes by unit of time
621    ! ------------------------------------------------
622
623    d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/dtime
624    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/dtime
625    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/dtime
626    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
627
628    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/dtime
629
630   print *,'zdu ', zdu
631   print *,'zdv ', zdv
632   print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1)
633
634    d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/dtime
635    d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/dtime
636    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/dtime
637    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/dtime
638
639    d_h_col = (zh_col(:,2)-zh_col(:,1))/dtime
640
641  end if ! end if (fl_ebil .GT. 0)
642!
643
644  RETURN
645END SUBROUTINE diag_phys_tend
646
647SUBROUTINE integr_v(nlon, nlev, zcpvap, &
648                    temp, qv, ql, qs, uu, vv, zairm,  &
649                    zqw_col, zql_col, zqs_col, zek_col, zh_dair_col, &
650                    zh_qw_col, zh_ql_col, zh_qs_col, zh_col)
651
652IMPLICIT none
653  include "YOMCST.h"
654
655INTEGER,                    INTENT(IN)    :: nlon,nlev
656REAL,                       INTENT(IN)    :: zcpvap
657REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: temp, qv, ql, qs, uu, vv
658REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: zairm
659REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqw_col
660REAL, DIMENSION(nlon),      INTENT(OUT)   :: zql_col
661REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqs_col
662REAL, DIMENSION(nlon),      INTENT(OUT)   :: zek_col
663REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_dair_col
664REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qw_col
665REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_ql_col
666REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qs_col
667REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_col
668
669INTEGER          :: i, k
670
671
672  ! Reset variables
673    zqw_col(:) = 0.
674    zql_col(:) = 0.
675    zqs_col(:) = 0.
676    zek_col(:) = 0.
677    zh_dair_col(:) = 0.
678    zh_qw_col(:) = 0.
679    zh_ql_col(:) = 0.
680    zh_qs_col(:) = 0.
681
682!JLD    write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
683 
684    ! ------------------------------------------------
685    ! Compute vertical sum for each atmospheric column
686    ! ------------------------------------------------
687    DO k = 1, nlev
688      DO i = 1, nlon
689        ! Watter mass
690        zqw_col(i) = zqw_col(i) + qv(i, k)*zairm(i, k)
691        zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
692        zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
693        ! Kinetic Energy
694        zek_col(i) = zek_col(i) + 0.5*(uu(i,k)**2+vv(i,k)**2)*zairm(i, k)
695        ! Air enthalpy : dry air, water vapour, liquid, solid
696        zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-qv(i,k)-ql(i,k)-qs(i,k))* &
697                                               zairm(i, k)*temp(i, k)
698        zh_qw_col(i) = zh_qw_col(i) +  zcpvap*temp(i, k)         *qv(i, k)*zairm(i, k)    !jyg
699        zh_ql_col(i) = zh_ql_col(i) + (zcpvap*temp(i, k) - rlvtt)*ql(i, k)*zairm(i, k)   !jyg
700        zh_qs_col(i) = zh_qs_col(i) + (zcpvap*temp(i, k) - rlstt)*qs(i, k)*zairm(i, k)   !jyg
701      END DO
702    END DO
703    ! compute total air enthalpy
704    zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:)
705
706END SUBROUTINE integr_v
707
708SUBROUTINE prt_enerbil (text, itap)
709!======================================================================
710! Print enenrgy budget diagnotics for the 1D case
711!======================================================================
712
713!======================================================================
714! Declarations
715!======================================================================
716
717USE dimphy, ONLY: klon, klev
718USE phys_state_var_mod, ONLY : dtime
719USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con
720USE geometry_mod, ONLY: longitude_deg, latitude_deg
721USE print_control_mod, ONLY: prt_level
722USE cmp_seri_mod
723USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
724  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
725USE phys_local_var_mod, ONLY: evap, sens
726USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri &
727   &  , rain_lsc, snow_lsc
728USE climb_hq_mod, ONLY : d_h_col_vdf, f_h_bnd
729IMPLICIT none
730include "YOMCST.h"
731
732! Arguments :
733!------------
734CHARACTER*(*) text ! text specifing the involved parametrization
735integer itap        ! time step number
736! local variables
737! ---------------
738real bilq_seuil,  bilh_seuil ! thresold on error in Q and H budget
739real bilq_error,  bilh_error ! erros in Q and H budget
740real bilq_bnd,  bilh_bnd     ! Q and H budget due to exchange with boundaries
741integer bilq_ok,  bilh_ok
742CHARACTER*(12) status
743
744bilq_seuil = 1.E-10
745bilh_seuil = 1.E-1
746bilq_ok=0
747bilh_ok=0
748
749!!print *,'prt_level:',prt_level,' fl_ebil:',fl_ebil,' fl_cor_ebil:',fl_cor_ebil
750if ( (fl_ebil .GT. 0) .and. (klon .EQ. 1)) then
751
752  bilq_bnd = 0.
753  bilh_bnd = 0.
754
755  param: SELECT CASE (text)
756  CASE("vdf") param
757      bilq_bnd = evap(1)
758      bilh_bnd = sens(1)+(rcpv-rcpd)*evap(1)*t_seri(1,1)
759  CASE("lsc") param
760      bilq_bnd = - rain_lsc(1) - snow_lsc(1)
761      bilh_bnd = (-(rcw-rcpd)*t_seri(1,1) + rlvtt) * rain_lsc(1) &
762    &         + (-(rcs-rcpd)*t_seri(1,1) + rlstt) * snow_lsc(1)
763  CASE("convection") param
764      bilq_bnd = - rain_con(1) - snow_con(1)
765      bilh_bnd = (-(rcw-rcpd)*t_seri(1,1) + rlvtt) * rain_con(1) &
766    &         + (-(rcs-rcpd)*t_seri(1,1) + rlstt) * snow_con(1)
767  CASE("SW") param
768      bilh_bnd = topsw(1) - solsw(1)
769  CASE("LW") param
770      bilh_bnd = -(toplw(1) + sollw(1))
771  CASE DEFAULT param
772      bilq_bnd = 0.
773      bilh_bnd = 0.
774  END SELECT param
775
776  bilq_error = d_qt_col(1) - bilq_bnd
777  bilh_error = d_h_col(1) - bilh_bnd
778! are the errors too large?
779  if ( abs(bilq_error) .gt. bilq_seuil) bilq_ok=1
780  if ( abs(bilh_error) .gt. bilh_seuil) bilh_ok=1
781!
782! Print diagnostics
783! =================
784  if ( (bilq_ok .eq. 0).and.(bilh_ok .eq. 0) ) then
785    status="enerbil-OK"
786  else
787    status="enerbil-PB"
788  end if
789
790  if ( prt_level .GE. 3) then
791    write(*,9010) text,status," itap:",itap,"enerbilERROR: Q", bilq_error,"  H", bilh_error
7929010  format (1x,A8,2x,A12,A6,I4,A18,E15.6,A5,E15.6)
793  end if
794  if ( prt_level .GE. 3) then
795    write(*,9000) text,"enerbil: Q,H,KE budget", d_qt_col(1), d_h_col(1),d_ek_col(1)
796  end if
797  if ( prt_level .GE. 5) then
798    write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd
799    write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1)
800    write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1)
801  end if
802
803  specific_diag: SELECT CASE (text)
804  CASE("vdf") specific_diag
805    if ( prt_level .GE. 5) then
806      write(*,9000) text,"enerbil: d_h, bilh, sens,t_seri", d_h_col(1), bilh_bnd, sens(1), t_seri(1,1)
807      write(*,9000) text,"enerbil: d_h_col_vdf, f_h, diff",d_h_col_vdf, f_h_bnd, bilh_bnd-sens(1)
808    end if
809  CASE("lsc") specific_diag
810    if ( prt_level .GE. 5) then
811      write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_lsc(1), rlvtt * rain_lsc(1), -(rcw-rcpd)*t_seri(1,1) * rain_lsc(1)
812      write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt * snow_lsc(1), -(rcs-rcpd)*t_seri(1,1) * snow_lsc(1)
813    end if
814  CASE("convection") specific_diag
815    if ( prt_level .GE. 5) then
816      write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt * rain_con(1), -(rcw-rcpd)*t_seri(1,1) * rain_con(1)
817      write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt * snow_con(1), -(rcs-rcpd)*t_seri(1,1) * snow_con(1)
818    end if
819  END SELECT specific_diag
820
8219000 format (1x,A8,2x,A35,10E15.6)
822
823end if ! end if (fl_ebil .GT. 0)
824
825END SUBROUTINE prt_enerbil
826
827END MODULE add_phys_tend_mod
Note: See TracBrowser for help on using the repository browser.