source: LMDZ5/branches/testing/libf/phylmd/add_phys_tend.F90 @ 2435

Last change on this file since 2435 was 2435, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2396:2434 into testing branch
  • 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: 8.7 KB
Line 
1!
2! $Id: add_phys_tend.F90 2435 2016-01-28 16:02:13Z fairhead $
3!
4SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text,abortphy)
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, ONLY: klon, klev
17USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, &
18                              t_seri
19USE phys_state_var_mod, ONLY: ftsol
20USE geometry_mod, ONLY: longitude_deg, latitude_deg
21USE print_control_mod, ONLY: prt_level
22IMPLICIT none
23  include "YOMCST.h"
24  include "clesphys.h"
25
26! Arguments :
27!------------
28REAL zdu(klon,klev),zdv(klon,klev)
29REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev),zdqi(klon,klev)
30REAL paprs(klon,klev+1)
31CHARACTER*(*) text
32INTEGER abortphy
33
34! Local :
35!--------
36REAL zt,zq
37REAL zq_int, zqp_int, zq_new
38
39REAL zqp(klev)
40
41INTEGER i, k,j
42INTEGER jadrs(klon*klev), jbad
43INTEGER jqadrs(klon*klev), jqbad
44INTEGER kadrs(klon*klev)
45INTEGER kqadrs(klon*klev)
46
47LOGICAL done(klon)
48
49integer debug_level
50logical, save :: first=.true.
51!$OMP THREADPRIVATE(first)
52INTEGER, SAVE :: itap
53!$OMP THREADPRIVATE(itap)
54!======================================================================
55! Initialisations
56
57      IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele
58                              ! a deja plante.
59
60     debug_level=10
61     if (first) then
62        itap=0
63        first=.false.
64     endif
65! Incrementer le compteur de la physique
66     itap   = itap + 1
67!======================================================================
68! Ajout des tendances sur le vent et l'eau liquide
69!======================================================================
70
71     u_seri(:,:)=u_seri(:,:)+zdu(:,:)
72     v_seri(:,:)=v_seri(:,:)+zdv(:,:)
73     ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
74     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
75
76!======================================================================
77! On ajoute les tendances de la temperature et de la vapeur d'eau
78! en verifiant que ca ne part pas dans les choux
79!======================================================================
80
81      jbad=0
82      jqbad=0
83      DO k = 1, klev
84         DO i = 1, klon
85            zt=t_seri(i,k)+zdt(i,k)
86            zq=q_seri(i,k)+zdq(i,k)
87            IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then
88            jbad = jbad + 1
89            jadrs(jbad) = i
90            kadrs(jbad) = k
91            ENDIF
92            IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
93            jqbad = jqbad + 1
94            jqadrs(jqbad) = i
95            kqadrs(jqbad) = k
96            ENDIF
97            t_seri(i,k)=zt
98            q_seri(i,k)=zq
99         ENDDO
100      ENDDO
101
102!=====================================================================================
103! Impression et stop en cas de probleme important
104!=====================================================================================
105
106IF (jbad .GT. 0) THEN
107      DO j = 1, jbad
108         i=jadrs(j)
109         if(prt_level.ge.debug_level) THEN
110          print*,'PLANTAGE POUR LE POINT i lon lat =',&
111                 i,longitude_deg(i),latitude_deg(i),text
112          print*,'l    T     dT       Q     dQ    '
113          DO k = 1, klev
114             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
115          ENDDO
116          call print_debug_phys(i,debug_level,text)
117         endif
118      ENDDO
119ENDIF
120!
121!=====================================================================================
122! Impression, warning et correction en cas de probleme moins important
123!=====================================================================================
124IF (jqbad .GT. 0) THEN
125      done(:) = .false.                         !jyg
126      DO j = 1, jqbad
127        i=jqadrs(j)
128          if(prt_level.ge.debug_level) THEN
129           print*,'WARNING  : EAU POUR LE POINT i lon lat =',&
130                  i,longitude_deg(i),latitude_deg(i),text
131           print*,'l    T     dT       Q     dQ    '
132           DO k = 1, klev
133              write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
134           ENDDO
135          endif
136          IF (ok_conserv_q) THEN
137!jyg<20140228 Corrections pour conservation de l'eau
138            IF (.NOT.done(i)) THEN                  !jyg
139              DO k = 1, klev
140                zqp(k) = max(q_seri(i,k),1.e-15)
141              ENDDO
142              zq_int  = 0.
143              zqp_int = 0.
144              DO k = 1, klev
145                zq_int  = zq_int  + q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
146                zqp_int = zqp_int + zqp(k)     *(paprs(i,k)-paprs(i,k+1))/Rg
147              ENDDO
148              if(prt_level.ge.debug_level) THEN
149               print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
150                                    i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
151              endif
152              DO k = 1, klev
153                zq_new = zqp(k)*zq_int/zqp_int
154                zdq(i,k) = zdq(i,k) + zq_new - q_seri(i,k)
155                q_seri(i,k) = zq_new
156              ENDDO
157              done(i) = .true.
158            ENDIF !(.NOT.done(i))
159          ELSE
160!jyg>
161            DO k = 1, klev
162              zq=q_seri(i,k)+zdq(i,k)
163              if (zq.lt.1.e-15) then
164                 if (q_seri(i,k).lt.1.e-15) then
165                  if(prt_level.ge.debug_level) THEN
166                   print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
167                  endif
168                  q_seri(i,k)=1.e-15
169                  zdq(i,k)=(1.e-15-q_seri(i,k))
170                 endif
171              endif
172!              zq=q_seri(i,k)+zdq(i,k)
173!              if (zq.lt.1.e-15) then
174!                 zdq(i,k)=(1.e-15-q_seri(i,k))
175!              endif
176            ENDDO
177!jyg<20140228
178          ENDIF   !  (ok_conserv_q)
179!jyg>
180      ENDDO ! j = 1, jqbad
181ENDIF
182!
183
184!IM ajout memes tests pour reverifier les jbad, jqbad beg
185      jbad=0
186      jqbad=0
187      DO k = 1, klev
188         DO i = 1, klon
189            IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
190            jbad = jbad + 1
191            jadrs(jbad) = i
192!            if(prt_level.ge.debug_level) THEN
193!             print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
194!            endif
195            ENDIF
196            IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
197            jqbad = jqbad + 1
198            jqadrs(jqbad) = i
199            kqadrs(jqbad) = k
200!            if(prt_level.ge.debug_level) THEN
201!             print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
202!            endif
203            ENDIF
204         ENDDO
205      ENDDO
206IF (jbad .GT. 0) THEN
207      DO j = 1, jbad
208         i=jadrs(j)
209         k=kadrs(j)
210         if(prt_level.ge.debug_level) THEN
211          print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',&
212                 i,itap,longitude_deg(i),latitude_deg(i),text,jbad, &
213       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
214!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
215          print*,'l    T     dT       Q     dQ    '
216          DO k = 1, klev
217             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
218          ENDDO
219          call print_debug_phys(i,debug_level,text)
220         endif
221      ENDDO
222ENDIF
223!
224IF (jqbad .GT. 0) THEN
225      DO j = 1, jqbad
226         i=jqadrs(j)
227         k=kqadrs(j)
228         if(prt_level.ge.debug_level) THEN
229          print*,'WARNING  : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',&
230                 i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,&
231       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
232!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
233          print*,'l    T     dT       Q     dQ    '
234          DO k = 1, klev
235            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
236          ENDDO
237          call print_debug_phys(i,debug_level,text)
238         endif
239      ENDDO
240ENDIF
241
242
243!======================================================================
244! Contrôle des min/max pour arrêt du modèle
245! Si le modele est en mode abortphy, on retire les tendances qu'on
246! vient d'ajouter. Pas exactement parce qu'on ne tient pas compte des
247! seuils.
248!======================================================================
249
250      CALL hgardfou(t_seri,ftsol,text,abortphy)
251      IF (abortphy==1) THEN
252        Print*,'ERROR ABORT hgardfou dans ',text
253        u_seri(:,:)=u_seri(:,:)-zdu(:,:)
254        v_seri(:,:)=v_seri(:,:)-zdv(:,:)
255        ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
256        qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
257      ENDIF
258
259
260
261      RETURN
262      END
Note: See TracBrowser for help on using the repository browser.