source: LMDZ6/trunk/libf/phylmd/add_phys_tend_mod.F90 @ 5134

Last change on this file since 5134 was 5051, checked in by evignon, 5 months ago

petit changement dans add_phys_tend suite au renommage des tendances de neige soufflee

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