source: LMDZ6/branches/Test_modipsl/libf/phylmd/add_phys_tend_mod.F90 @ 5201

Last change on this file since 5201 was 4523, checked in by evignon, 19 months ago

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