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

Last change on this file since 2791 was 2641, checked in by Laurent Fairhead, 8 years ago

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