source: LMDZ5/trunk/libf/phylmd/add_phys_tend_mod.F90 @ 2836

Last change on this file since 2836 was 2812, checked in by jyg, 8 years ago

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