source: LMDZ5/trunk/libf/phylmd/add_phys_tend.F90 @ 2545

Last change on this file since 2545 was 2400, checked in by fhourdin, 9 years ago

Compilation bug fixing ...
  • 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
RevLine 
[904]1!
[1163]2! $Id: add_phys_tend.F90 2400 2015-11-21 18:38:32Z fairhead $
[904]3!
[2235]4SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text,abortphy)
[904]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
[2399]16USE dimphy, ONLY: klon, klev
17USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, &
18                              t_seri
[2400]19USE phys_state_var_mod, ONLY: ftsol
[2399]20USE geometry_mod, ONLY: longitude_deg, latitude_deg
21USE print_control_mod, ONLY: prt_level
[904]22IMPLICIT none
[2007]23  include "YOMCST.h"
24  include "clesphys.h"
[904]25
26! Arguments :
27!------------
28REAL zdu(klon,klev),zdv(klon,klev)
[2086]29REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev),zdqi(klon,klev)
[1998]30REAL paprs(klon,klev+1)
[904]31CHARACTER*(*) text
[2235]32INTEGER abortphy
[904]33
34! Local :
35!--------
36REAL zt,zq
[1998]37REAL zq_int, zqp_int, zq_new
[904]38
[1998]39REAL zqp(klev)
40
[904]41INTEGER i, k,j
[924]42INTEGER jadrs(klon*klev), jbad
43INTEGER jqadrs(klon*klev), jqbad
[972]44INTEGER kadrs(klon*klev)
45INTEGER kqadrs(klon*klev)
[904]46
[1998]47LOGICAL done(klon)
48
[904]49integer debug_level
[972]50logical, save :: first=.true.
[987]51!$OMP THREADPRIVATE(first)
[972]52INTEGER, SAVE :: itap
[987]53!$OMP THREADPRIVATE(itap)
[904]54!======================================================================
55! Initialisations
56
[2235]57      IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele
58                              ! a deja plante.
59
60     debug_level=10
[972]61     if (first) then
62        itap=0
63        first=.false.
64     endif
65! Incrementer le compteur de la physique
66     itap   = itap + 1
[904]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(:,:)
[2086]74     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
[904]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
[972]90            kadrs(jbad) = k
[904]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
[972]95            kqadrs(jqbad) = k
[904]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)
[1047]109         if(prt_level.ge.debug_level) THEN
[2399]110          print*,'PLANTAGE POUR LE POINT i lon lat =',&
111                 i,longitude_deg(i),latitude_deg(i),text
[1047]112          print*,'l    T     dT       Q     dQ    '
113          DO k = 1, klev
[904]114             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
[1047]115          ENDDO
116          call print_debug_phys(i,debug_level,text)
117         endif
[904]118      ENDDO
119ENDIF
120!
121!=====================================================================================
122! Impression, warning et correction en cas de probleme moins important
123!=====================================================================================
124IF (jqbad .GT. 0) THEN
[1998]125      done(:) = .false.                         !jyg
[904]126      DO j = 1, jqbad
[1998]127        i=jqadrs(j)
128          if(prt_level.ge.debug_level) THEN
[2399]129           print*,'WARNING  : EAU POUR LE POINT i lon lat =',&
130                  i,longitude_deg(i),latitude_deg(i),text
[1998]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
[2007]136          IF (ok_conserv_q) THEN
[1998]137!jyg<20140228 Corrections pour conservation de l'eau
[2007]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))
[1998]175!              endif
[2007]176            ENDDO
177!jyg<20140228
178          ENDIF   !  (ok_conserv_q)
[1998]179!jyg>
[2007]180      ENDDO ! j = 1, jqbad
[904]181ENDIF
182!
183
[972]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
[1163]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
[972]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
[1163]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
[972]203            ENDIF
204         ENDDO
205      ENDDO
206IF (jbad .GT. 0) THEN
207      DO j = 1, jbad
208         i=jadrs(j)
209         k=kadrs(j)
[1047]210         if(prt_level.ge.debug_level) THEN
[2399]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, &
[972]213       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
[1047]214!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
[972]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
[1163]222ENDIF
[972]223!
224IF (jqbad .GT. 0) THEN
225      DO j = 1, jqbad
226         i=jqadrs(j)
227         k=kqadrs(j)
[1047]228         if(prt_level.ge.debug_level) THEN
[2399]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,&
[972]231       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
[1047]232!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
[972]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
[2235]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
[904]261      RETURN
262      END
Note: See TracBrowser for help on using the repository browser.