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

Last change on this file since 4384 was 4384, checked in by Sebastien Nguyen, 17 months ago

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