source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/add_phys_tend.F90 @ 3817

Last change on this file since 3817 was 3817, checked in by millour, 10 years ago

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