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

Last change on this file since 2128 was 2086, checked in by fhourdin, 10 years ago

nclusion de la thermodynamique de la glace
Ice thermodynamics
(Catherine Rio)
  • 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: 7.7 KB
RevLine 
[904]1!
[1163]2! $Id: add_phys_tend.F90 2086 2014-07-09 21:19:07Z lguez $
[904]3!
[2086]4SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text)
[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
16use dimphy
17use phys_local_var_mod
[941]18use phys_state_var_mod
[904]19IMPLICIT none
[2007]20  include "iniprint.h"
21  include "YOMCST.h"
22  include "clesphys.h"
[904]23
24! Arguments :
25!------------
26REAL zdu(klon,klev),zdv(klon,klev)
[2086]27REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev),zdqi(klon,klev)
[1998]28REAL paprs(klon,klev+1)
[904]29CHARACTER*(*) text
30
31! Local :
32!--------
33REAL zt,zq
[1998]34REAL zq_int, zqp_int, zq_new
[904]35
[1998]36REAL zqp(klev)
37
[904]38INTEGER i, k,j
[924]39INTEGER jadrs(klon*klev), jbad
40INTEGER jqadrs(klon*klev), jqbad
[972]41INTEGER kadrs(klon*klev)
42INTEGER kqadrs(klon*klev)
[904]43
[1998]44LOGICAL done(klon)
45
[904]46integer debug_level
[972]47logical, save :: first=.true.
[987]48!$OMP THREADPRIVATE(first)
[972]49INTEGER, SAVE :: itap
[987]50!$OMP THREADPRIVATE(itap)
[904]51!======================================================================
52! Initialisations
53
54debug_level=10
[972]55     if (first) then
56        itap=0
57        first=.false.
58     endif
59! Incrementer le compteur de la physique
60     itap   = itap + 1
[904]61!======================================================================
62! Ajout des tendances sur le vent et l'eau liquide
63!======================================================================
64
65     u_seri(:,:)=u_seri(:,:)+zdu(:,:)
66     v_seri(:,:)=v_seri(:,:)+zdv(:,:)
67     ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
[2086]68     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
[904]69
70!======================================================================
71! On ajoute les tendances de la temperature et de la vapeur d'eau
72! en verifiant que ca ne part pas dans les choux
73!======================================================================
74
75      jbad=0
76      jqbad=0
77      DO k = 1, klev
78         DO i = 1, klon
79            zt=t_seri(i,k)+zdt(i,k)
80            zq=q_seri(i,k)+zdq(i,k)
81            IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then
82            jbad = jbad + 1
83            jadrs(jbad) = i
[972]84            kadrs(jbad) = k
[904]85            ENDIF
86            IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
87            jqbad = jqbad + 1
88            jqadrs(jqbad) = i
[972]89            kqadrs(jqbad) = k
[904]90            ENDIF
91            t_seri(i,k)=zt
92            q_seri(i,k)=zq
93         ENDDO
94      ENDDO
95
96!=====================================================================================
97! Impression et stop en cas de probleme important
98!=====================================================================================
99
100IF (jbad .GT. 0) THEN
101      DO j = 1, jbad
102         i=jadrs(j)
[1047]103         if(prt_level.ge.debug_level) THEN
104          print*,'PLANTAGE POUR LE POINT i rlon rlat =',i,rlon(i),rlat(i),text
105          print*,'l    T     dT       Q     dQ    '
106          DO k = 1, klev
[904]107             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
[1047]108          ENDDO
109          call print_debug_phys(i,debug_level,text)
110         endif
[904]111      ENDDO
112ENDIF
113!
114!=====================================================================================
115! Impression, warning et correction en cas de probleme moins important
116!=====================================================================================
117IF (jqbad .GT. 0) THEN
[1998]118      done(:) = .false.                         !jyg
[904]119      DO j = 1, jqbad
[1998]120        i=jqadrs(j)
121          if(prt_level.ge.debug_level) THEN
122           print*,'WARNING  : EAU POUR LE POINT i rlon rlat =',i,rlon(i),rlat(i),text
123           print*,'l    T     dT       Q     dQ    '
124           DO k = 1, klev
125              write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
126           ENDDO
127          endif
[2007]128          IF (ok_conserv_q) THEN
[1998]129!jyg<20140228 Corrections pour conservation de l'eau
[2007]130            IF (.NOT.done(i)) THEN                  !jyg
131              DO k = 1, klev
132                zqp(k) = max(q_seri(i,k),1.e-15)
133              ENDDO
134              zq_int  = 0.
135              zqp_int = 0.
136              DO k = 1, klev
137                zq_int  = zq_int  + q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/Rg
138                zqp_int = zqp_int + zqp(k)     *(paprs(i,k)-paprs(i,k+1))/Rg
139              ENDDO
140              if(prt_level.ge.debug_level) THEN
141               print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
142                                    i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
143              endif
144              DO k = 1, klev
145                zq_new = zqp(k)*zq_int/zqp_int
146                zdq(i,k) = zdq(i,k) + zq_new - q_seri(i,k)
147                q_seri(i,k) = zq_new
148              ENDDO
149              done(i) = .true.
150            ENDIF !(.NOT.done(i))
151          ELSE
152!jyg>
153            DO k = 1, klev
154              zq=q_seri(i,k)+zdq(i,k)
155              if (zq.lt.1.e-15) then
156                 if (q_seri(i,k).lt.1.e-15) then
157                  if(prt_level.ge.debug_level) THEN
158                   print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
159                  endif
160                  q_seri(i,k)=1.e-15
161                  zdq(i,k)=(1.e-15-q_seri(i,k))
162                 endif
163              endif
164!              zq=q_seri(i,k)+zdq(i,k)
165!              if (zq.lt.1.e-15) then
166!                 zdq(i,k)=(1.e-15-q_seri(i,k))
[1998]167!              endif
[2007]168            ENDDO
169!jyg<20140228
170          ENDIF   !  (ok_conserv_q)
[1998]171!jyg>
[2007]172      ENDDO ! j = 1, jqbad
[904]173ENDIF
174!
175
[972]176!IM ajout memes tests pour reverifier les jbad, jqbad beg
177      jbad=0
178      jqbad=0
179      DO k = 1, klev
180         DO i = 1, klon
181            IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
182            jbad = jbad + 1
183            jadrs(jbad) = i
[1163]184!            if(prt_level.ge.debug_level) THEN
185!             print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
186!            endif
[972]187            ENDIF
188            IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
189            jqbad = jqbad + 1
190            jqadrs(jqbad) = i
191            kqadrs(jqbad) = k
[1163]192!            if(prt_level.ge.debug_level) THEN
193!             print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
194!            endif
[972]195            ENDIF
196         ENDDO
197      ENDDO
198IF (jbad .GT. 0) THEN
199      DO j = 1, jbad
200         i=jadrs(j)
201         k=kadrs(j)
[1047]202         if(prt_level.ge.debug_level) THEN
203          print*,'PLANTAGE2 POUR LE POINT i itap rlon rlat txt jbad zdt t',i,itap,rlon(i),rlat(i),text,jbad, &
[972]204       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
[1047]205!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
[972]206          print*,'l    T     dT       Q     dQ    '
207          DO k = 1, klev
208             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
209          ENDDO
210          call print_debug_phys(i,debug_level,text)
211         endif
212      ENDDO
[1163]213ENDIF
[972]214!
215IF (jqbad .GT. 0) THEN
216      DO j = 1, jqbad
217         i=jqadrs(j)
218         k=kqadrs(j)
[1047]219         if(prt_level.ge.debug_level) THEN
220          print*,'WARNING  : EAU2 POUR LE POINT i itap rlon rlat txt jqbad zdq q zdql ql',i,itap,rlon(i),rlat(i),text,jqbad,&
[972]221       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
[1047]222!!!       if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
[972]223          print*,'l    T     dT       Q     dQ    '
224          DO k = 1, klev
225            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
226          ENDDO
227          call print_debug_phys(i,debug_level,text)
228         endif
229      ENDDO
230ENDIF
231
232      CALL hgardfou(t_seri,ftsol,text)
[904]233      RETURN
234      END
Note: See TracBrowser for help on using the repository browser.