source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/add_phys_tend.F90 @ 5416

Last change on this file since 5416 was 2007, checked in by fhourdin, 11 years ago

Modification pour garantir la conservation des espèces traces.

Modification for tracer conservation

Jean-Yves Grandpeix
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.6 KB
Line 
1!
2! $Id: add_phys_tend.F90 2007 2014-04-06 12:20:38Z fairhead $
3!
4SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,paprs,text)
5!======================================================================
6! Ajoute les tendances des variables physiques aux variables
7! d'etat de la dynamique t_seri, q_seri ...
8! On en profite pour faire des tests sur les tendances en question.
9!======================================================================
10
11
12!======================================================================
13! Declarations
14!======================================================================
15
16use dimphy
17use phys_local_var_mod
18use phys_state_var_mod
19IMPLICIT none
20  include "iniprint.h"
21  include "YOMCST.h"
22  include "clesphys.h"
23
24! Arguments :
25!------------
26REAL zdu(klon,klev),zdv(klon,klev)
27REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev)
28REAL paprs(klon,klev+1)
29CHARACTER*(*) text
30
31! Local :
32!--------
33REAL zt,zq
34REAL zq_int, zqp_int, zq_new
35
36REAL zqp(klev)
37
38INTEGER i, k,j
39INTEGER jadrs(klon*klev), jbad
40INTEGER jqadrs(klon*klev), jqbad
41INTEGER kadrs(klon*klev)
42INTEGER kqadrs(klon*klev)
43
44LOGICAL done(klon)
45
46integer debug_level
47logical, save :: first=.true.
48!$OMP THREADPRIVATE(first)
49INTEGER, SAVE :: itap
50!$OMP THREADPRIVATE(itap)
51!======================================================================
52! Initialisations
53
54debug_level=10
55     if (first) then
56        itap=0
57        first=.false.
58     endif
59! Incrementer le compteur de la physique
60     itap   = itap + 1
61!======================================================================
62! Ajout des tendances sur le vent et l'eau liquide
63!======================================================================
64
65     u_seri(:,:)=u_seri(:,:)+zdu(:,:)
66     v_seri(:,:)=v_seri(:,:)+zdv(:,:)
67     ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
68
69!======================================================================
70! On ajoute les tendances de la temperature et de la vapeur d'eau
71! en verifiant que ca ne part pas dans les choux
72!======================================================================
73
74      jbad=0
75      jqbad=0
76      DO k = 1, klev
77         DO i = 1, klon
78            zt=t_seri(i,k)+zdt(i,k)
79            zq=q_seri(i,k)+zdq(i,k)
80            IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then
81            jbad = jbad + 1
82            jadrs(jbad) = i
83            kadrs(jbad) = k
84            ENDIF
85            IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
86            jqbad = jqbad + 1
87            jqadrs(jqbad) = i
88            kqadrs(jqbad) = k
89            ENDIF
90            t_seri(i,k)=zt
91            q_seri(i,k)=zq
92         ENDDO
93      ENDDO
94
95!=====================================================================================
96! Impression et stop en cas de probleme important
97!=====================================================================================
98
99IF (jbad .GT. 0) THEN
100      DO j = 1, jbad
101         i=jadrs(j)
102         if(prt_level.ge.debug_level) THEN
103          print*,'PLANTAGE POUR LE POINT i rlon rlat =',i,rlon(i),rlat(i),text
104          print*,'l    T     dT       Q     dQ    '
105          DO k = 1, klev
106             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
107          ENDDO
108          call print_debug_phys(i,debug_level,text)
109         endif
110      ENDDO
111ENDIF
112!
113!=====================================================================================
114! Impression, warning et correction en cas de probleme moins important
115!=====================================================================================
116IF (jqbad .GT. 0) THEN
117      done(:) = .false.                         !jyg
118      DO j = 1, jqbad
119        i=jqadrs(j)
120          if(prt_level.ge.debug_level) THEN
121           print*,'WARNING  : EAU POUR LE POINT i rlon rlat =',i,rlon(i),rlat(i),text
122           print*,'l    T     dT       Q     dQ    '
123           DO k = 1, klev
124              write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
125           ENDDO
126          endif
127          IF (ok_conserv_q) THEN
128!jyg<20140228 Corrections pour conservation de l'eau
129            IF (.NOT.done(i)) THEN                  !jyg
130              DO k = 1, klev
131                zqp(k) = max(q_seri(i,k),1.e-15)
132              ENDDO
133              zq_int  = 0.
134              zqp_int = 0.
135              DO k = 1, klev
136                zq_int  = zq_int  + q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
137                zqp_int = zqp_int + zqp(k)     *(paprs(i,k)-paprs(i,k+1))/Rg
138              ENDDO
139              if(prt_level.ge.debug_level) THEN
140               print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
141                                    i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
142              endif
143              DO k = 1, klev
144                zq_new = zqp(k)*zq_int/zqp_int
145                zdq(i,k) = zdq(i,k) + zq_new - q_seri(i,k)
146                q_seri(i,k) = zq_new
147              ENDDO
148              done(i) = .true.
149            ENDIF !(.NOT.done(i))
150          ELSE
151!jyg>
152            DO k = 1, klev
153              zq=q_seri(i,k)+zdq(i,k)
154              if (zq.lt.1.e-15) then
155                 if (q_seri(i,k).lt.1.e-15) then
156                  if(prt_level.ge.debug_level) THEN
157                   print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
158                  endif
159                  q_seri(i,k)=1.e-15
160                  zdq(i,k)=(1.e-15-q_seri(i,k))
161                 endif
162              endif
163!              zq=q_seri(i,k)+zdq(i,k)
164!              if (zq.lt.1.e-15) then
165!                 zdq(i,k)=(1.e-15-q_seri(i,k))
166!              endif
167            ENDDO
168!jyg<20140228
169          ENDIF   !  (ok_conserv_q)
170!jyg>
171      ENDDO ! j = 1, jqbad
172ENDIF
173!
174
175!IM ajout memes tests pour reverifier les jbad, jqbad beg
176      jbad=0
177      jqbad=0
178      DO k = 1, klev
179         DO i = 1, klon
180            IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
181            jbad = jbad + 1
182            jadrs(jbad) = i
183!            if(prt_level.ge.debug_level) THEN
184!             print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
185!            endif
186            ENDIF
187            IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
188            jqbad = jqbad + 1
189            jqadrs(jqbad) = i
190            kqadrs(jqbad) = k
191!            if(prt_level.ge.debug_level) THEN
192!             print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
193!            endif
194            ENDIF
195         ENDDO
196      ENDDO
197IF (jbad .GT. 0) THEN
198      DO j = 1, jbad
199         i=jadrs(j)
200         k=kadrs(j)
201         if(prt_level.ge.debug_level) THEN
202          print*,'PLANTAGE2 POUR LE POINT i itap rlon rlat txt jbad zdt t',i,itap,rlon(i),rlat(i),text,jbad, &
203       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
204!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
205          print*,'l    T     dT       Q     dQ    '
206          DO k = 1, klev
207             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
208          ENDDO
209          call print_debug_phys(i,debug_level,text)
210         endif
211      ENDDO
212ENDIF
213!
214IF (jqbad .GT. 0) THEN
215      DO j = 1, jqbad
216         i=jqadrs(j)
217         k=kqadrs(j)
218         if(prt_level.ge.debug_level) THEN
219          print*,'WARNING  : EAU2 POUR LE POINT i itap rlon rlat txt jqbad zdq q zdql ql',i,itap,rlon(i),rlat(i),text,jqbad,&
220       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
221!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
222          print*,'l    T     dT       Q     dQ    '
223          DO k = 1, klev
224            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
225          ENDDO
226          call print_debug_phys(i,debug_level,text)
227         endif
228      ENDDO
229ENDIF
230
231      CALL hgardfou(t_seri,ftsol,text)
232      RETURN
233      END
Note: See TracBrowser for help on using the repository browser.