source: LMDZ6/branches/Amaury_dev/libf/phylmd/add_phys_tend_mod.F90 @ 5136

Last change on this file since 5136 was 5134, checked in by abarral, 8 weeks ago

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