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

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

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