source: LMDZ6/trunk/libf/dyn3dmem/qminimum_loc.F @ 3877

Last change on this file since 3877 was 3801, checked in by lguez, 4 years ago

Bug fix from revision 3800
  • 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:keywords set to Id
File size: 10.0 KB
RevLine 
[3800]1!
2!     $Id: qminimum_loc.F 3801 2021-01-16 18:55:09Z acozic $
3!
[2270]4      SUBROUTINE qminimum_loc( q,nqtot,deltap )
[1823]5      USE parallel_lmdz
[3800]6      USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif,             &
7     &   ratiomin,qperemin ! CRisi 23nov2020
[1632]8      IMPLICIT none
9c
10c  -- Objet : Traiter les valeurs trop petites (meme negatives)
11c             pour l'eau vapeur et l'eau liquide
12c
[2600]13      include "dimensions.h"
14      include "paramet.h"
[3801]15      include "iniprint.h"
[1632]16c
[2270]17      INTEGER nqtot ! CRisi: on remplace nq par nqtot
18      REAL q(ijb_u:ije_u,llm,nqtot), deltap(ijb_u:ije_u,llm)
[1632]19c
20      INTEGER iq_vap, iq_liq
21      PARAMETER ( iq_vap = 1 ) ! indice pour l'eau vapeur
22      PARAMETER ( iq_liq = 2 ) ! indice pour l'eau liquide
23      REAL seuil_vap, seuil_liq
24      PARAMETER ( seuil_vap = 1.0e-10 ) ! seuil pour l'eau vapeur
25      PARAMETER ( seuil_liq = 1.0e-11 ) ! seuil pour l'eau liquide
26c
27c  NB. ....( Il est souhaitable mais non obligatoire que les valeurs des
28c            parametres seuil_vap, seuil_liq soient pareilles a celles
29c            qui  sont utilisees dans la routine    ADDFI       )
30c     .................................................................
31c
32      INTEGER i, k, iq
33      REAL zx_defau, zx_abc, zx_pump(ijb_u:ije_u), pompe
[2270]34
35      real zx_defau_diag(ijb_u:ije_u,llm,2)
36      real q_follow(ijb_u:ije_u,llm,2)
[1632]37c
38      REAL SSUM
39      EXTERNAL SSUM
40c
41      INTEGER imprim
42      SAVE imprim
43      DATA imprim /0/
44c$OMP THREADPRIVATE(imprim)
45      INTEGER ijb,ije
46      INTEGER Index_pump(ij_end-ij_begin+1)
47      INTEGER nb_pump
[2270]48      INTEGER ixt
49      INTEGER iso_verif_noNaN_nostop
[1632]50c
51c Quand l'eau liquide est trop petite (ou negative), on prend
52c l'eau vapeur de la meme couche et la convertit en eau liquide
53c (sans changer la temperature !)
54c
55
[3800]56        !write(lunout,*) 'qminimum 52: entree'
[2270]57        if (ok_iso_verif) then
58           call check_isotopes(q,ij_begin,ij_end,'qminimum 52')   
59        endif !if (ok_iso_verif) then     
60
[1632]61      ijb=ij_begin
62      ije=ij_end
63
[2270]64      zx_defau_diag(ijb:ije,:,:)=0.0
[2285]65      q_follow(ijb:ije,:,1:2)=q(ijb:ije,:,1:2) 
[2270]66
[3800]67      !write(lunout,*) 'qminimum 57'
[1632]68c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
69      DO 1000 k = 1, llm
70      DO 1040 i = ijb, ije
71            if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
[2270]72
73              if (ok_isotopes) then
74                 zx_defau_diag(i,k,iq_liq)=AMAX1
75     :               ( seuil_liq - q(i,k,iq_liq), 0.0 )
76              endif !if (ok_isotopes) then
77
[1632]78               q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
79               q(i,k,iq_liq) = seuil_liq
80            endif
81 1040 CONTINUE
82 1000 CONTINUE
83c$OMP END DO NOWAIT
84c$OMP BARRIER
85c --->  SYNCHRO OPENMP ICI
86
[2270]87
[1632]88c
89c Quand l'eau vapeur est trop faible (ou negative), on complete
90c le defaut en prennant de l'eau vapeur de la couche au-dessous.
91c
[3800]92      !write(lunout,*) 'qminimum 81'
[1632]93      iq = iq_vap
94c
95      DO k = llm, 2, -1
96ccc      zx_abc = dpres(k) / dpres(k-1)
97c$OMP DO SCHEDULE(STATIC)
98      DO i = ijb, ije
[2270]99
[1632]100         if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
[2270]101
102            if (ok_isotopes) then
103              zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 )
104            endif !if (ok_isotopes) then
105
[1632]106            q(i,k-1,iq) =  q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) *
107     &           deltap(i,k) / deltap(i,k-1)
108            q(i,k,iq)   =  seuil_vap 
[2270]109
[1632]110         endif
111      ENDDO
112c$OMP END DO NOWAIT
113      ENDDO
114c$OMP BARRIER
[2270]115
[1632]116c
117c Quand il s'agit de la premiere couche au-dessus du sol, on
118c doit imprimer un message d'avertissement (saturation possible).
119c
[3800]120      !write(lunout,*) 'qminimum 106'
[1632]121      nb_pump=0
122c$OMP DO SCHEDULE(STATIC)
123      DO i = ijb, ije
124         zx_pump(i) = AMAX1( 0.0, seuil_vap - q(i,1,iq) )
125         q(i,1,iq)  = AMAX1( q(i,1,iq), seuil_vap )
126         IF (zx_pump(i) > 0.0) THEN
127            nb_pump = nb_pump+1
128            Index_pump(nb_pump)=i
129         ENDIF
130      ENDDO
131c$OMP END DO 
132!      pompe = SSUM(ije-ijb+1,zx_pump(ijb),1)
133
134      IF (imprim.LE.100 .AND. nb_pump .GT. 0 ) THEN
135         PRINT *, 'ATT!:on pompe de l eau au sol'
136         DO i = 1, nb_pump
137               imprim = imprim + 1
138               PRINT*,'  en ',index_pump(i),zx_pump(index_pump(i))
139         ENDDO
140      ENDIF
[2270]141
[3800]142      !write(lunout,*) 'qminimum 128'
[2270]143      if (ok_isotopes) then
[3800]144              !write(lunout,*) 'qminimum 140'
[2270]145      ! CRisi: traiter de même les traceurs d'eau
146      ! Mais il faut les prendre à l'envers pour essayer de conserver la
147      ! masse.
148      ! 1) pompage dans le sol 
149      ! On suppose que ce pompage se fait sans isotopes -> on ne modifie
150      ! rien ici et on croise les doigts pour que ça ne soit pas trop
151      ! génant
[3800]152      ! en fait, si, c'est genant quand les isotopes doivent eux même transporter des
153      ! traceurs -> apporter aussi un peu d'isotopes... Combien?
154      ! Essayer tnat/2 = -500 permil? C'est déjà mieux que -1000
155      ! permil...
156      ! pb: que faire pour les traceurs?
157c$OMP DO SCHEDULE(STATIC)     
[2270]158      DO i = ijb, ije
159        if (zx_pump(i).gt.0.0) then
160          q_follow(i,1,iq_vap)=q_follow(i,1,iq_vap)+zx_pump(i)
161        endif !if (zx_pump(i).gt.0.0) then
162      enddo !DO i = ijb, ije 
[3800]163c$OMP END DO
[2270]164
165      ! 2) transfert de vap vers les couches plus hautes
[3800]166      !write(lunout,*) 'qminimum 158'
[2270]167      do k=2,llm
[3800]168c$OMP DO SCHEDULE(STATIC)     
[2270]169        DO i = ijb, ije
170          if (zx_defau_diag(i,k,iq_vap).gt.0.0) then             
[3800]171              ! on ajoute la vapeur en k     
172!              write(lunout,*) 'i,k,q_follow(i,k-1,iq_vap)=',
173!     :                 i,k,q_follow(i,k-1,iq_vap)         
174              if (q_follow(i,k-1,iq_vap).lt.qperemin) then
175                write(lunout,*) 'tmp qmin: on stoppe'
176                write(lunout,*) 'zx_pump(i)=',zx_pump(i)
177                write(lunout,*) 'q_follow(i,:,iq_vap)=',
178     :                   q_follow(i,:,iq_vap)
179                write(lunout,*) 'k=',k
180                call abort_gcm("qminimum","not enough vapor",1)
181              endif 
182            do ixt=1,ntraciso
183!                write(lunout,*) 'qmin 168: ixt=',ixt
184!                write(lunout,*) 'q(i,k,iqiso(ixt,iq_vap)=',
185!     :             q(i,k,iqiso(ixt,iq_vap))
186!                write(lunout,*) 'zx_defau_diag(i,k,iq_vap)=',
187!     :                  zx_defau_diag(i,k,iq_vap)
188!                write(lunout,*) 'q(i,k-1,iqiso(ixt,iq_vap)=',
189!     :                   q(i,k-1,iqiso(ixt,iq_vap))     
190
[2270]191               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
192     :              +zx_defau_diag(i,k,iq_vap)
193     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
194               
195              if (ok_iso_verif) then
196                if (iso_verif_noNaN_nostop(q(i,k,iqiso(ixt,iq_vap)),
197     :                   'qminimum 155').eq.1) then
198                   write(*,*) 'i,k,ixt=',i,k,ixt
199                   write(*,*) 'q_follow(i,k-1,iq_vap)=',
200     :                   q_follow(i,k-1,iq_vap)
201                   write(*,*) 'q(i,k,iqiso(ixt,iq_vap))=',
202     :                   q(i,k,iqiso(ixt,iq_vap))
203                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
204     :                   zx_defau_diag(i,k,iq_vap)
205                   write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=',
206     :                   q(i,k-1,iqiso(ixt,iq_vap))
207                   stop
208                endif
209              endif
210
211              ! et on la retranche en k-1
212               q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap))
213     :              -zx_defau_diag(i,k,iq_vap)
214     :              *deltap(i,k)/deltap(i,k-1)
215     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
216
217               if (ok_iso_verif) then
218                if (iso_verif_noNaN_nostop(q(i,k-1,iqiso(ixt,iq_vap)),
219     :                   'qminimum 175').eq.1) then
220                   write(*,*) 'k,i,ixt=',k,i,ixt
221                   write(*,*) 'q_follow(i,k-1,iq_vap)=',
222     :                   q_follow(i,k-1,iq_vap)
223                   write(*,*) 'q(i,k,iqiso(ixt,iq_vap))=',
224     :                   q(i,k,iqiso(ixt,iq_vap))
225                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
226     :                   zx_defau_diag(i,k,iq_vap)
227                   write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=',
228     :                   q(i,k-1,iqiso(ixt,iq_vap))
229                   stop
230                endif
231              endif
232
233              enddo !do ixt=1,niso
234              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
235     :               +zx_defau_diag(i,k,iq_vap)
236              q_follow(i,k-1,iq_vap)=   q_follow(i,k-1,iq_vap)
237     :               -zx_defau_diag(i,k,iq_vap)
238     :              *deltap(i,k)/deltap(i,k-1)
239          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
240        enddo !DO i = 1, ip1jmp1       
[3800]241c$OMP END DO
242        enddo !do k=2,llm
[2270]243
244        if (ok_iso_verif) then
245           call check_isotopes(q,ijb,ije,'qminimum 168')
246        endif !if (ok_iso_verif) then
247       
248     
249        ! 3) transfert d'eau de la vapeur au liquide
250        !write(*,*) 'qminimum 164'
251        do k=1,llm
[3800]252c$OMP DO SCHEDULE(STATIC)
[2270]253        DO i = ijb, ije
254          if (zx_defau_diag(i,k,iq_liq).gt.0.0) then
255
256              ! on ajoute eau liquide en k en k             
257              do ixt=1,ntraciso
258               q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq))
259     :              +zx_defau_diag(i,k,iq_liq)
260     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)
261              ! et on la retranche à la vapeur en k
262               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
263     :              -zx_defau_diag(i,k,iq_liq)
264     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)   
265              enddo !do ixt=1,niso
266              q_follow(i,k,iq_liq)=   q_follow(i,k,iq_liq)
267     :               +zx_defau_diag(i,k,iq_liq)
268              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
269     :               -zx_defau_diag(i,k,iq_liq)
270          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
[3800]271        enddo !DO i = ijb, ije
272c$OMP END DO       
[2270]273       enddo !do k=2,llm 
274
275        if (ok_iso_verif) then
276           call check_isotopes(q,ijb,ije,'qminimum 197')
277        endif !if (ok_iso_verif) then
278
279      endif !if (ok_isotopes) then
280      !write(*,*) 'qminimum 188'
[1632]281c
282      RETURN
283      END
Note: See TracBrowser for help on using the repository browser.