source: LMDZ6/trunk/libf/phylmdiso/add_phys_tend_mod.F90 @ 5419

Last change on this file since 5419 was 5285, checked in by abarral, 8 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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