source: LMDZ6/trunk/libf/phylmd/add_phys_tend_mod.F90 @ 3639

Last change on this file since 3639 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

  • as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

  • due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

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