source: LMDZ6/trunk/libf/phylmd/add_phys_tend_mod.f90 @ 5451

Last change on this file since 5451 was 5285, checked in by abarral, 2 months ago

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