source: LMDZ6/trunk/libf/phylmdiso/reevap.F90 @ 4009

Last change on this file since 4009 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 7.0 KB
RevLine 
[3927]1  SUBROUTINE reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2   &         d_t_eva,d_q_eva,d_ql_eva,d_qs_eva &
3#ifdef ISO
4             ,xt_seri,xtl_seri,xts_seri,d_xt_eva,d_xtl_eva,d_xts_eva &
5#endif
6   &     )
7
8    ! flag to include modifications to ensure energy conservation (if flag >0)
9    USE add_phys_tend_mod, only : fl_cor_ebil
10#ifdef ISO
11    USE infotrac_phy, ONLY: ntraciso   
12#ifdef ISOVERIF
13    USE isotopes_verif_mod
14!, ONLY: errmax,errmaxrel, iso_verif_o18_aberrant_nostop,deltaD,deltaO
15    USE isotopes_mod, ONLY: iso_eau,iso_hdo,iso_o18,ridicule
16#ifdef ISOTRAC
17    USE isotrac_routines_mod, ONLY: iso_verif_traceur_pbidouille   
18#endif
19#endif
20#endif
21    IMPLICIT none
22    !>======================================================================
23
24    INTEGER klon,klev,iflag_ice_thermo
25    REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri,q_seri,ql_seri,qs_seri
26    REAL, DIMENSION(klon,klev), INTENT(out) :: d_t_eva,d_q_eva,d_ql_eva,d_qs_eva
27
28    REAL za,zb,zdelta,zlvdcp,zlsdcp
29    INTEGER i,k
30
31#ifdef ISO
32    REAL, DIMENSION(ntraciso,klon,klev), INTENT(in) :: xt_seri,xtl_seri,xts_seri
33    REAL, DIMENSION(ntraciso,klon,klev), INTENT(out) :: d_xt_eva,d_xtl_eva,d_xts_eva
34    integer ixt
35#endif
36
37    !--------Stochastic Boundary Layer Triggering: ALE_BL--------
38    !---Propri\'et\'es du thermiques au LCL
39    include "YOMCST.h"
40    include "YOETHF.h"
41    include "FCTTRE.h"
42    !IM 100106 BEG : pouvoir sortir les ctes de la physique
43    !
44    ! Re-evaporer l'eau liquide nuageuse
45    !
46!print *,'rrevap ; fl_cor_ebil:',fl_cor_ebil,' iflag_ice_thermo:',iflag_ice_thermo,' RVTMP2',RVTMP2
47    DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
48       DO i = 1, klon
49        if (fl_cor_ebil .GT. 0) then
50          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
51          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
52        else
53          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
54          !jyg<
55          !  Attention : Arnaud a propose des formules completement differentes
56          !                  A verifier !!!
57          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
58        end if
59          IF (iflag_ice_thermo .EQ. 0) THEN
60             zlsdcp=zlvdcp
61          ENDIF
62          !>jyg
63
64          IF (iflag_ice_thermo.eq.0) THEN   
65             !pas necessaire a priori
66
67             zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
68  zdelta = 0.
69             zb = MAX(0.0,ql_seri(i,k))
70             za = - MAX(0.0,ql_seri(i,k)) &
71                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
72             d_t_eva(i,k) = za
73             d_q_eva(i,k) = zb
74             d_ql_eva(i,k) = -ql_seri(i,k)
75             d_qs_eva(i,k) = 0.
76
77#ifdef ISO
78         do ixt=1,ntraciso
79            zb = MAX(0.0,xtl_seri(ixt,i,k))
80            d_xt_eva(ixt,i,k) = zb
81            d_xtl_eva(ixt,i,k) = -xtl_seri(ixt,i,k)
82            d_xts_eva(ixt,i,k) = 0.
83         enddo ! do ixt=1,ntraciso
84#ifdef ISOVERIF
85      do ixt=1,ntraciso
86        call iso_verif_noNaN(xt_seri(ixt,i,k), &
87     &     'physiq 2417: apres evap tot')
88      enddo
89      if (iso_eau.gt.0) then
90              call iso_verif_egalite_choix( &
91     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
92     &          'physiq 1891+, après reevap totale',errmax,errmaxrel)
93              call iso_verif_egalite_choix( &
94     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
95     &          'physiq 2209+, après reevap totale',errmax,errmaxrel)
96       endif !if (iso_eau.gt.0) then       
97      if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
98            if (q_seri(i,k).gt.ridicule) then 
99               if (iso_verif_o18_aberrant_nostop( &
100     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
101     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
102     &              'physiq 2315: apres reevap totale').eq.1) then
103                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)
104                  write(*,*) 'd_q_eva(i,k)=',d_q_eva(i,k)
105                  write(*,*) 'deltaD(d_q_eva(i,k))=',deltaD(d_xt_eva(iso_HDO,i,k)/d_q_eva(i,k))
106                  write(*,*) 'deltaO18(d_q_eva(i,k))=',deltaO(d_xt_eva(iso_O18,i,k)/d_q_eva(i,k))
107                  stop
108              endif !  if (iso_verif_o18_aberrant_nostop
109            endif !if (q_seri(i,k).gt.errmax) then   
110        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
111#ifdef ISOTRAC     
112             call iso_verif_traceur(xt_seri(1,i,k), &
113     &           'physiq 2165')
114             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
115     &           'physiq 2165b')
116#endif               
117         
118#endif           
119#endif
120
121          ELSE
122
123             !CR: on r\'e-\'evapore eau liquide et glace
124
125             !        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
126             !        zb = MAX(0.0,ql_seri(i,k))
127             !        za = - MAX(0.0,ql_seri(i,k)) &
128             !             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
129             zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
130             za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
131                  - MAX(0.0,qs_seri(i,k))*zlsdcp
132             d_t_eva(i,k) = za
133             d_q_eva(i,k) = zb
134             d_ql_eva(i,k) = -ql_seri(i,k)
135             d_qs_eva(i,k) = -qs_seri(i,k)
136
137#ifdef ISO
138         do ixt=1,ntraciso
139            zb = MAX(0.0,xtl_seri(ixt,i,k)+xts_seri(ixt,i,k))
140            d_xt_eva(ixt,i,k) = zb
141            d_xtl_eva(ixt,i,k) = -xtl_seri(ixt,i,k)
142            d_xts_eva(ixt,i,k) = -xts_seri(ixt,i,k)
143         enddo ! do ixt=1,ntraciso
144
145#ifdef ISOVERIF
146      do ixt=1,ntraciso
147      call iso_verif_noNaN(xt_seri(ixt,i,k), &
148     &     'physiq 2417: apres evap tot')
149      enddo
150      if (iso_eau.gt.0) then
151              call iso_verif_egalite_choix( &
152     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
153     &          'physiq 1891, après réévap totale',errmax,errmaxrel)
154              call iso_verif_egalite_choix( &
155     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
156     &          'physiq 2209, après réévap totale',errmax,errmaxrel)
157              call iso_verif_egalite_choix( &
158     &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
159     &          'physiq 2209b, après réévap totale',errmax,errmaxrel)
160       endif !if (iso_eau.gt.0) then
161     
162      if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
163            if (q_seri(i,k).gt.ridicule) then 
164               if (iso_verif_o18_aberrant_nostop( &
165     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
166     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
167     &              'physiq 2408: apres reevap totale').eq.1) then
168                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)                       
169                  stop
170              endif !  if (iso_verif_o18_aberrant_nostop
171            endif !if (q_seri(i,k).gt.errmax) then   
172        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then       
173#ifdef ISOTRAC     
174             call iso_verif_traceur(xt_seri(1,i,k), &
175     &           'physiq 2165')
176             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
177     &           'physiq 2165b')
178#endif                 
179#endif                 
180#endif
181
182          ENDIF
183
184       ENDDO
185    ENDDO
186
187RETURN
188
189END SUBROUTINE reevap
Note: See TracBrowser for help on using the repository browser.