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

Last change on this file since 4077 was 4004, checked in by Laurent Fairhead, 3 years ago

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